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SUMMARY 


The manner of arriving at the stability determinant for nat 
ural motions of a bending vs. torsion aeroelastic system is out- 
This permits determination of the characteristics of the 
Sample 


lined. 
natural modes of the system at any chosen air speed. 
calculations are presented for a series of wings subject to bending 
vs. torsion flutter, showing how the natural mode dampings and 
frequencies change as the flutter speed is approached. It is 
found that the damping in the critical mode increases with air 
speed until at least 85 per cent of the critical speed is attained; it 
then rapidly decreases to the flutter value of zero. The non 
critical mode encounters a sharp increase in stability in the prox 
imity of flutter. A comparison is also made between the normal 
mode characteristics and the system response due to forced ex 
citation; it is found that the latter usually affords an earlier 
indication of the approach to a critical state. From this it ap 
pears that flight-flutter testing is best carried out using the forced 


vibration technique. 


INTRODUCTION 


— ei—maramead TREND in current aeronautical engi- 
neering practice is the growing detail with which 


a variety of dynamical problems must be treated. 
Whereas in the past it has been possible to approach 
a design almost solely on the basis of satisfying certain 
pseudostatic operating conditions, when modern high- 
performance aircraft are involved, instances 
arise where this point of view must be discarded; in its 
place, a realistic consideration of the craft as a free 
This 
trend is particularly evident in certain problems in the 
fields of stability and structures, where the complexity 
of questions relating to handling qualities, maneuver- 
ability, and dynamic loads has increased to the point 


many 


body in three-dimensional motion is necessary. 
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where the simplified pseudostatic approximations no 
longer apply. Examples of such problems are found in 
the study of the aircraft behavior in gusty air, the 
dynamical loads induced during landing, and the ef- 
fects of aeroelasticity during high-performance opera- 
tion. 

The present paper is concerned with the dynamic 
characteristics of an aircraft while in flight at an air 
speed close to a critical wing bending-torsion flutter 
speed. The study is of interest for several reasons: 
The nature of the dynamic behavior of the craft near 
flutter is of importance in deciding on the maximum 
allowable flight speed consistent with safe operation. 
Also, one procedure for experimentally determining the 
flutter speed of a craft is by the method of forced os- 
cillation;! to utilize this procedure with safety, the 
performance near flutter should be thoroughly under- 
Finally, the results obtained here shed light on 
the mechanism of the bending-torsion 
flutter mode. Since it appears that many high-per- 
formance aircraft will of necessity operate at speeds 
that are relatively close toa flutter speed, questions of 
the kind discussed here assume an increased significance. 

Considering now an elastic aircraft in free flight, it 
follows that the motions of the craft result from the 
simultaneous action of inertia, elastic, and aerodynamic 
forces, all of which are a consequence of the motion of 
the craft itself; if only small departures from a pseudo- 
static state of motion are considered, as is the usual 
case in theoretical aeroelastic treatments, then the 
magnitudes of these forces are linearly dependent on 
From the general theory 


stood. 
aeroelastic 


the disturbances themselves. 
for linear dynamic systems, the result then appears that 
the disturbed motion is compounded from a number of 
natural modes, whose descriptions follow from a knowl- 
edge of the roots of the pertinent stability determinant. 
In addition to the forces mentioned, dissipative forces 
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are encountered in all physical systems; while these are 
of considerable practical importance, their influence on 
the dynamic behavior of the system is well understood, 
and it is permissible to neglect their effects for the 
purposes of the present theoretical discussion. 

The nature of the response mechanism for the case of 
forced oscillations in flight is also disclosed by the gen- 
eral theory. The response is composed of two parts: 
one is the natural response arising from the comple- 
mentary solution to the equations of motion and in- 
volving the natural modes, while the second part is 
the forced response and corresponds to the particular 
solution to the equations. 

The present paper is devoted to a study of the natural 
modes only. The character of the forced response has 
received the attention of numerous previous authors,” * 
of which the cited references are examples. Myklestad?* 
derives the forced response resulting from the applica- 
tion of an external force of harmonically oscillating mag- 
nitude, while Smilg and Wasserman® investigate the 
hypothetical structural damping required to maintain 
constant amplitude oscillations at chosen air speeds. 
Studies of this kind, while informative in their own 
right, do not disclose the character of the natural modes. 

In the present work, the general stability determi 
nant for an aeroelastic system involving wing bending 
vs. wing torsion is established. The equations of 
motion are derived by making use of the Wagner air 
forces,* © which presume a potential flow of an in- 
compressible fluid past a thin, two-dimensional airfoil 
in unsteady motion. It is then shown how the roots of 
the stability determinant can be extracted in a rela- 
tively simple manner, enabling a description of the 
natural modes. Sample calculations for a number of 
aircraft wings are presented over speed ranges from 
zero air speed up to the flutter speed. Finally, the re- 
sults of the natural mode study are compared with the 
character of the forced response under conditions of 
harmonic oscillation to determine to what extent the 
two display a common mechanism. 

The preparation of the present paper was suggested 
by certain findings by the authors in connection with a 
program sponsored by the Dynamics Branch, Air- 
craft Laboratory, Air Materiel Command, U.S. Air 
Force. Messrs. Smilg, Spielberg, and Wasserman, all 
of that agency, have contributed toward the develop- 


ment of many of the ideas. 


DERIVATION OF THE AEROELASTIC MODES 


Since the concern here is with the natural modes at 
air speeds below the flutter speed, it follows that the 
analysis must be of sufficient generality to embrace 
nonharmonic motions of the system; for example, it is 
to be expected that the natural modes will be character- 
ized by decaying oscillatory motions, and the air forces 
introduced in the equations of motion must embrace 


this circumstance. While structural damping is not 
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included in the work here for reasons already given, it is 
of passing interest to note that no physically substan- 
tiated description of the dissipation mechanism for gen- 
eral kinds of motion is at present available to the 
analyst. 

The Wagner air forces,* *® which apply under the 
conditions mentioned in the previous section of the 
paper, are employed in the present analysis. For the 
convenience of the reader, a brief description of these 
air forces is now given. Throughout the analysis in 
this paper, a nondimensional coordinate system is used, 
according to the following scheme: Physical lengths 
are referred to the reference length 5,, where 0, is a 
reference semichord length of the wing; physical times 
are referred to b,/Uo, where Up is the forward speed of 
the craft, considered to remain constant in the present 
study; physical masses are referred to pb,*, where p is 
the mass density of the air; physical forces are then 
referred to the force pl/y’b,’. 

Consider now an airfoil with a semichord 6 moving 
horizontally to the left with a uniform forward velocity 
in a two-dimensional field of flow. Let w,(7) be the time 
history of the downward velocity of the airfoil meas- 
ured at the rear quarter-chord point, and let @(r) be 
the time history of the nose-up pitching motions (r is 
the nondimensional time): The forces and moments 
acting per unit of span as a result of unsteady motion 
are then as follows: 

(a) A/ift force with center of pressure at the forward 


quarter-chord point and of amount 


Li = ono (7 . *) ta(71) + O(71)} dr; + 
0 ) 


7 
2arbk () }w, (0) + 6(0); ( 
) 


where the function «(7) describes the growth of lift on 
an airfoil of physical semichord 6, due to a sudden (step 
function) change in angle of attack and is equal to 


K(r) = 1 — ®(7) 


where ®(7) is the Wagner ‘‘lift deficiency”’ function cal- 
culated by von Karman and Sears.® Jones® suggests 
that x(r) can be approximated with sufficient accuracy 
for practical purposes by the expression 


Bar (33) 


k(r) = 1 — ae °"” — ace” 


where a; = 0.165, a2 = 0.335, B, = 0.0455, and B: = 
0.3. For convenience, this approximation is retained 
in the work reported here. 

(b) A /ift force with center of pressure at the mid- 


chord and of amount 
Le = mb*{w,(7r) + 6(7r)} (4 


(c) A downward force with center of pressure hallf- 
way between the forward quarter-chord and the mid- 
chord and of amount 


L3 = (2/2)b3{6(r)} (0 
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(d) Finally, a mose-down moment of amount 


M, = (w/2)b3{6(r)} (6) 


It is evident that L; describes the circulatory compo- 
nent of lift; Zz and Z; are “‘apparent mass” contribu- 
tions, while 4/7, is a rotational damping component. 

The total lift on the airfoil per unit of span is then 


iF = Li + Le —_ Ls 


(7) 
and the total moment about the forward quarter-chord, 
positive when nose-up, is 


M = —(1/2)L. + (1/4)L3; — M, (3) 


In the analysis to follow, interest is centered on the 
air forces that arise from a damped, oscillatory motion 
of the airfoil. A general motion of this kind is defined 
by the time histories 

s= Ae’; @ = Be’; (9) 


A= -a+ jw 


where 3 is a complex coordinate, the real (or imaginary) 
part of which represents the vertical displacement at 
the forward quarter-chord point; 6 is now taken to be 
a complex coordinate, the real (or imaginary) part of 
which describes the pitch of the airfoil, positive when 
nose-up; and A and B are complex constants. The 
parameter a describes the degree of the decay of the 
motion, while w refers to the circular frequency; evi- 
dently, by proper choice of a and w, the time histories 
can be made to apply to harmonic, as well as to aperi- 
odic, motions. It follows at once that w,(r) can be 
represented in complex form as 


w,(r) = 5 + 66 = A(z + DO) (10) 


The air forces due to the motion Eq. (9) are readily 
obtained from an application of Eqs. (1) through (10). 
Treating first the circulatory component, this is most 
easily derived with the aid of operational methods. 
The Laplace transform of L; is shown by the standard 


procedures to be? 
6p" + bop + 63 
(bp + B:)(bp + Be)(p — AX) 
[XA + (b+ 1)B] 


Ll, = 2xb 
(11) 


where p is the Laplace operator, 


= [1 — o, — o)d 
bo = [8i(1 — ds) + Bo(1 a 4) |b 
63 = BiBo 


For the problem at hand, the aperiodic air force com- 
ponents arising from the poles p = —£,/b and p = 
— 62/6 are not of interest; only the air forces due to the 
pole = \ are required, and these are 


(d:A? + der + 43) 


Li = 2x5 —— Az + (6A + 196] (12) 
(br + By) (bX - Bo) ) | : 

Itis convenient to set 
c(X) (6,A? oa boA + 63) [(br oa Bi) (br + Bo) | (13) 
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and then the circulatory component of the wing lift 
becomes 


L, = 2nb{c(d)} [Az + (OA + 16] (14) 


The remaining lift and moment components Le, Ls, W/; 
are derived without difficulty by direct substitution of 
the zs, 6 forms. 

It is then found that the total wing lift and the wing 
moment about the forward quarter-chord point, each 
per unit of span, can be written as 


L=12+ 10 (15) 
M = m,z + m6 (16) 

where 
1, = wb[2dc(d) + ddr?*] (17) 
ly = wb[2(bX + 1)c(A) + OA(DA + 1) — (6?/2)A*] (18) 
m, = —(m/2)b*d? (19) 


e 1 b r 
Me = rh? | — 2 (DA + 1) + 54 ; — 1] (20) 


Before concluding this résumé of the air forces in- 
duced by the unsteady motion of the airfoil, it is of in- 
terest to note how the present descriptions reduce to the 
well-known results for certain limiting cases. For the 
pseudostatic case, when \ = 0, the familiar expressions 
for the air forces arise (note that the static angle of 
attack is Az + 0). Considering next harmonic motion 
with the circular frequency Q, so that a = 0 and \ = 
jQ, it is seen that Eqs. (15) through (20) reduce to the 
Theodorsen results,* * with the function c(A) fulfilling 
the role of his C(k) function. If Eq. (13) is examined 
with this in mind, it is found that the function c(d) 
agrees to within a few per cent with the corresponding 
values of C(k) over the entire range of values of 2. 

Attention is now given to the derivation of the equa- 
tions of motion for the bending vs. torsion aeroelastic 
system. Consider first a two-dimensional system of 
unit span. Letting z, be the downward displacement 
at the elastic axis, then the kinetic energy for the system 
can be written as 


2T = mz? + 283.0 + 16? 


(21) 


where m is the airfoil mass, S is the mass static unbal- 
ance about the elastic axis (positive for a rearward 
c.g.), and J is the mass moment of inertia about the 
elastic axis. The potential energy of the system is 
written as 

2V = ma,2.2 + [w,?0? (22) 
where w, and w, are the uncoupled wing bending and 
torsion frequencies, respectively. The air forces, Eqs. 
(15) and (16), are next referenced to the elastic axis 
displacement, so that the lift and the nose-up moment 
about the elastic axis are written as 


i = 12, f 1,0 
M, = + m6 


(23) 


mS, (24) 
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Fic. 3. 
the wing bending-torsion system with ws 
and « = 304 cycles per min 





Natural mode characteristics plotted vs. 


Natural mode characteristics plotted vs. 





air speed for 
= 103.5 cycles per min 
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where 
l; = i 1, = le — el, 2.) 
m, = m. + el, (26) 
Mz = My + e(ly — m.) — e7l, (27) 


and where e is the rearward distance from the quarter- 
chord point to the elastic axis. 

An application of the Lagrange equations of motion, 
using Eqs. (9) and (21) through (24) then leads to the 
equations of motion for the system, 


}m(A2 + wy?) + Ate, + (SA2 +4386 =O (28) 


1}SA2 — mts. + L(A? + w,2) — mo}@ = 0 (29) 


Those values of \ which are the roots of the stability 
determinant arising from Eqs. (28) and (29) describe 
the natural modes of the aeroelastic system. With the 
roots \ known, the mode shapes are determined from 


the equations of motion in the usual manner. 


In the 
far been confined to a two-dimensional system. 


consideration has thus 
The 


interests of clarity, 


however, are derived for a wing strip of 


air forces, 
arbitrary chord 6, so that the tools necessary for a 
three-dimensional aeroelastic analysis are at hand. 


The procedure for extending the argument to a three- 
dimensional system is well-known® and is not repeated 
here. The sample calculations described later are car- 
ried out for a three-dimensional wing. 


METHOD OF SOLUTION 


For the purposes of the present study, attention is 
centered on those aeroelastic modes that have their 
zero air-speed origins in the two ground vibration 
modes that are composed of fundamental wing bending 
and torsion. As the air speed is gradually increased, 
the character of these modes also changes in a continu- 
ous fashion, so that a space plot of A vs. Up can be vis 
ualized as having two branches: The ‘“‘bending branch” 
has its zero air-speed origin in the ground mode which 
is primarily wing bending, while the “torsion branch” 
has an analogous interpretation. 

Since the parameters /,, /2, 77, and m2 are themselves 
complex functions of \, it is evident that obtaining a 
straightforward solution for \ from the stability de- 
terminant of Eqs. (28) and (29) is a tedious operation; 
it is necessary that the roots of an eighth degree poly- 
nomial be found. Fortunately, however, the amount 
labor required to effect a solution is considerably 
an iterative procedure for evalu- 


reduced by employing 
This procedure is described 


ating the stability roots. 
as follows: 

It is to be noted that the numerical values of w, and 
w, are dependent on the magnitude of Up; 
separate solution must be carried out for each air speed. 
To start a particular solution, attention is centered on 
one of the modes, say the bending branch, and a rea- 
sonable estimate for the mode frequency is made. 
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air speed for the wing bending-torsion system with w, = 


500 





and 


207 cycles per min. 


w = 304 cycles per min. 


fe) 
big. 4. Natural mode characteristics plotted vs. 
Assuming that a = 0 and setting w equal to the esti- 


mated frequency [see Eq. (9)], the associated value of 
c(\) is computed from Eq. (13), and then /;, /2, mm, and 
m, are evaluated from Eqs. (25) through (27). The 
trial values for the latter four coefficients are next intro- 
duced into the stability determinant, which now takes 
the form of a quadratic in \*. The roots of the quad- 
ratic found, and a refined value for the 

sponding to the bending branch is thus obtained; 


corre- 
four 
roots \ are obtained from the quadratic, and the root 


are 


of interest is readily chosen by inspection and on the 
basis of physical reasoning. 

Based on the refined \ value, a new value for c(A) is 
and then the four air-force coefficients /;, 
The new quadratic in \° 


computed, 
lo, my, Wt2 are recalculated. 
corresponding to these coefficients is written and the 
roots derived. A second refinement for the value of the 
bending branch A is thus obtained. The process is re- 
peated until the values of } computed by successive 
agree with each other within reasonable 
i.e., the process is continued until convergence 


iterations 
limits 
is attained. In order to compute the torsion branch 
\, the entire iteration process is repeated from the be- 
ginning. 

It is found that convergence is usually attained within 
a surprisingly few iteration cycles, three or four such 
cycles generally being sufficient for each \ determina- 
tion. The numerical operations involved in the work 
are relatively convenient to handle with desk-type 
calculating machines, so that solutions can be obtained 
rapidly and without tedious effort. 


CALCULATIONS FOR SAMPLE WINGS 


To illustrate the kind of results to be expected from 
an aeroelastic study of the type described here, a series 
ol calculations are now presented which bring to light 





the effect of a change in the ratio w/w, on the aero 
elastic behavior of a fighter-class wing. A series of 
wings are studied, each being dynamically identical with 
the next except that the natural frequency of the nor 
mal bending mode is varied. The frequency of the 
normal (uncoupled) torsion mode is kept at 304 cycles 
and wings with the following (uncoupled) 
dealt with: 69, 103.5, 
Bi corresponding 
227, 0.341, 0.454, 


).224, 
The wings are idiom as supported 


per min., 
bending mode frequencies are 
138, 207, and 276 cycles per min. 
ratios of w,/w, for the wings 
0.682, and 0.909. 
in the “‘free-free”’ 
and a substantial concentrated weight at the wing tip 
(arising from a tip tank installation). The detailed 
data for the wing system with the 138 cycles per min. 


are 


sense and carry a fuselage at the root 


bending frequency are drawn from the description of 
an actual model of a modern fighter aircraft, the system 
being fully described in a previous paper by the pres- 
ent authors;® see the Case I aircraft, Table I, of that 
All of the wings dealt with here are assumed 
to have the same generalized masses as the 138 cycle 


paper. 


per min. wing, as well as the same geometric and mode 
Hence, the wings are dynami- 
for the 


shape characteristics. 


cally identical, save bending frequency pa- 


rameter. 

The results of the calculations for the bending branch 
and torsion branch stability roots for the five wings 
are shown in Figs. 1-5. The calculations extend from 
zero air speed up to the flutter speed, and the plots 
The 
air speeds quoted are uncorrected for compressibility 
effects—i.e., they refer to flight in 
air. 


are of @ and w as a function of the air speed Up. 
‘“incompressible”’ 


The first point of interest to be noted in Figs. 1—5 con- 
cerns the variation of the mode damping a@ with air 
For the branch that leads to flutter, it is seen 


speed. 
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Fic. 5. Natural mode characteristics plotted vs. air speed for 
the wing bending-torsion system with w, = 276 cycles per min. 
and w = 304 cycles per min 


that in each of the five cases the mode damping at first 
gradually increases with air speed; however, at some 
point in the region of from 85 to 95 per cent of the flutter 
speed, the mode stability suddenly starts to decrease 
and the damping rapidly drops to the flutter value 
(a = 0). For the mode that does not lead to flutter, 
the damping gradually increases with the air speed until 
the region near flutter is reached. Then, as the critical 
speed is approached and the one mode goes to flutter, 
the damping in the second mode suddenly increases 
at a rapid rate. This substantial improvement in 
the stability of the nonflutter mode as the critical 
speed is approached is a matter of considerable inter- 
est. 

Experience with wing bending-torsion flutter has 
shown it to be a precipitous phenomenon, the transition 
from stability to instability occurring within a narrow 
speed range; it is to be noted that this type of behavior 
is indicated by the present calculations. Also, it is a 
generally accepted design specification within the in- 
dustry that the maximum flight speed for a model shall 
not exceed 85 per cent of its flutter speed. For wing 
bending-torsion flutter this is seen to be a reasonable 
requirement, since a decrease in mode stability is not to 
be expected until this limit speed is surpassed. 

Turning now to a consideration of the mode fre- 
quencies, it follows from Figs. 1-5 that the bending 
branch frequency remains approximately the same and 
equal to the ground frequency at all air speeds, although 
a sudden slight increase is observed in the immediate 
vicinity of the flutter speed. The frequency of the 
torsion branch mode falls appreciably, however, as the 
air speed increases; at the flutter speed, the frequencies 
of both branches lie fairly close to each other for all five 
cases investigated. As in the three other curves in 
each figure, a sudden change in the nature of the tor- 
sion branch frequency variation with air speed occurs 
when the flutter speed is approached. 

It should be noted that, for the systems with a low 
w»,/ w, ratio, it is the torsion branch mode that leads to 
flutter; for the higher ratios of the bending to torsion 
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frequencies, the mode connected with the bending 
branch results in flutter. The well-known drop in 
critical speed as the bending and torsion frequencies 
approach each other in value is also observed in the 
present calculations. 


res 


As a final point of interest, the logarithmic decre- 
ments for the critical modes in the five cases are com 
pared. It is readily shown that the logarithmic de- 
crement for the natural oscillations at a particular air 


speed is given by 
log A = 2r(a/w) 


Compare now the logarithmic decrements for the 
critical modes for the five wings at that air speed where 
the value of a is a maximum; this corresponds closely 
to the air speed at which a maximum of mode stability is 
indicated. In the order of ascending bending fre- 
quency, the values are roughly 1.08, 1.06, 0.67, 0.21 

and 0.074. 


It is thus seen that, as the ratio w,/w, increases, the 


















rate of decay of the natural oscillations in the critical , 
mode tends to decrease—i.e., the mode gives the ap- ‘ 
pearance of relatively lesser stability. As w/w, ap- 
proaches unity, the natural oscillations disappear at ‘ 
an extremely slow rate, regardless of air speed. Hence, le 
the severity of flutter decreases as the normal bend- d 
ing and torsion frequencies approach each other in (c 
value. pr 
st 

THE BEHAVIOR NEAR FLUTTER: A COMPARISON WITH in 
THE AIR MATERIEL COMMAND AND MYKLESTAD ca 
METHODS th 

A proper analysis of the behavior of an aeroelastic " 
system in the immediate vicinity of flutter is of con- es 
siderable importance when, for example, a program of “4 

flight-flutter testing is undertaken. The results of the 

analysis must be used to appraise the safety of the pro- I 
cedure in general and must also serve as a basis for ” 
fixing a safe increment of air speed for successive = 
tests.' In the simplest kind of flight-flutter testing, si 
the aircraft is disturbed from steady flight by the pilot, ol 
and the resulting rate of decay is observed; here, only du 
the natural modes are of interest. For more elaborate for 
tests employing forced vibration, the nature of the = 
forced response is of primary interest, although the crit 
presence of a poorly damped natural mode adds to the 1 
difficulty of obtaining accurate test data. forc 
It is now of interest to compare the behavior of the usu 
natural modes near flutter with the responses that re fut 
sult from forced excitation, the primary object being to the 
determine which of the two mechanisms is best adapted sale 
to flight-flutter studies. As has already been mentioned, : 
the Air Materiel Command method?’ investigates oath 
the damping required to sustain steady oscillations at date 
the various air speeds, while Myklestad? studies the of t 
response induced by an oscillatory applied force. The data 
results afforded by these two procedures are now com for ¢ 
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Fic. 6. Curves of g vs. Up calculated by the A.M.C. method 
for the wing bending-torsion systems with w, = 69, 138, and 207 
cycles per min. and a = 304 cycles per min 


pared with those obtained from an analysis of the na- 
tural modes as outlined in this paper. 

Fig. 6 shows the Uy vs. g curves derived by the A.M.C. 
method for three of the five wings studied in the last 
section (g is the fictitious damping coefficient required 
for steady oscillation). Only the curve for the critical 
mode in each case is shown. Comparing first the curve 
in Fig. 6 for the «, = 69 cycle per min. wing with the 
critical a vs. Up curve from Fig. 1, it is seen that the 
latter indicates a much sharper onset of flutter than 
does the former. For the w, = 138 cycle per min. wing 
(see Figs. 3 and 6), a similar situation prevails; the ap- 
proach to flutter is more severe in the natural mode 
study than in the Uy vs. g curve. For this case it is 
interesting to note that the natural mode analysis indi- 
cates that the bending branch leads to flutter, whereas 
the A.M.C. method shows the torsion branch becom- 
ing unstable. Finally, for the w, = 207 cycle per min. 
wing (Figs. 4 and 6), the two analyses predict about 
equal severity in the flutter approach.* 

Fig. 7 is taken from a flutter analysis of the w, = 
138 cycles per min. wing, carried out by the Myklestad 
It shows the imaginary part of the general- 
ized complex bending force required to sustain steady 


method. 


oscillations at various air speeds and at several values 
of the reduced speed; flutter is encountered at a re- 
duced speed of 16.75. It is seen that the shaking 
force starts to drop toward the flutter value (zero) 
at a somewhat earlier relative air speed than does the 
critical mode damping data shown in Fig. 6. 

These limited results tend to show that the use of the 
forced excitation method for flight-flutter testing will 
usually afford an earlier indication of the approach to 
Hence, 
the former method is to be preferred on the grounds of 


flutter than will a study of the natural modes. 
safety. 
* It is to be noted that the flutter speeds indicated in Figs. 1-5 


The 
data for the last two figures were calculated using the exact values 


do not coincide exactly with those shown in Figs. 6 and 7. 


of the flutter aerodynamic coefficients, whereas the Figs. 1-5 
data are subject to the errors arising from the approximation used 
for ¢(X), see Eq. (13). 
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CONCLUSIONS 


The present paper outlines the manner of arriving at 
the stability determinant for an aeroelastic system 
which permits the determination of the natural mode 
characteristics at any chosen flight speed. Although a 
wing bending vs. torsion system is specifically dealt 
with, extension of the analysis to other arrangements 
presents no theoretic difficulties. 

Sample calculations for a series of wing systems sub 
ject to bending-torsion flutter show that the damping in 
the critical natural mode increases until at least 85 per 
cent of the flutter speed is attained; it then rapidly 
drops to the flutter value of zero as the critical air speed 
isapproached. The maximum damping encountered in 
the critical mode decreases as the ratio of the normal 
bending and torsion frequencies in the system ap 
proaches unity. 

For the noncritical mode, a sharp increase in damping 
is indicated as the air speed comes in close proximity to 
the critical speed. For systems with a ratio of normal 
bending to torsion frequencies less than one, the fre- 
quency of the bending branch mode does not change 
markedly as flutter is approached; however, that of the 
torsion branch mode drops so that both modes have 
frequencies close to each other. 

Comparing the behavior of the natural modes near 
flutter with the system response induced by forced 
excitation, it appears that the forced response usually 
gives the earliest indication of an approach to a critical 
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Fic. 7. Complex part of shaking force required to sustain 
steady oscillations, plotted vs. air speed and for several values of 
reduced speed, as calculated by the Myklestad method for the 
wing bending-torsion system with w, = 138 cycles per min. and 
«, = 304 cycles per min. (Flutter is at a reduced velocity of 


16.75.) 





396 OF THE 


JOURNAL 
state. Hence, for flight-flutter testing of bending- 
torsion configurations, the method of forced excitation 


usually affords greater safety than a study of the na- 
tural modes. 
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The Flutter of Servo-Controlled Aircraft 


JONATHAN WINSON* 


Fairchild Engine and Airplane Corporation 


SUMMARY 


This paper considers the problem of the flutter of servo-con- 
trolled aircraft. The existence of flutter in such craft depends, 
in part, on the dynamic characteristics of the servomechanisms 
used for control. In classical flutter theory, the control surface is 
assumed to be part of a simple elastic system made up of inertia, 
spring, and damper components. With this representation, the 
equations of motion for control surface-wing motion are written. 
Their solution yields the flutter information required. In the 
servo-controlled aircraft, control surface motion under so-called 
locked stick conditions is an expression of the reaction of the 
highly complex servomechanism to external disturbance. The 
mechanism is usually nonlinear and difficult to treat mathemat- 
ically. 

In the present solution for the flutter of the servo-controlled 
craft, experimentally determined characteristics of the servo unit 
are used. A method is evolved wherein the impedance of the 
servo and related elastic structure, as obtained by test, is com- 
pared to a set of theoretical flutter boundary curves. The flutter 
boundary curves are functions of the aircraft parameters. The 
comparison yields a complete understanding of the flutter char- 
acteristics of the aircraft. The technique involved is general, 
embracing classical flutter as a special case. It is shown that the 
method is a revealing means of solving the classical flutter prob- 
lem. A sample calculation is carried out for a particular ship 
with and without servo-control. 


INTRODUCTION 


ier mameirgatead OF VARIOUS TYPES have achieved 
widespread use in the control of aircraft. Their 
use is universal in the intrinsically high-speed pilotless 
planes. With the development of servo-control has 
come the need to re-examine the concept of the flutter 
of an aircraft member in light of the contribution of the 
servomechanism to the flutter system. 

As applied to conventional aircraft, flutter is a self- 
excited instability of a physical system involving a 
flexible wing (with or without a control surface) and the 
air through which the wing travels. To make possible 
the mathematical treatment of this system, it is as- 
sumed that the motion of the wing is compounded of a 
number of modes, each distinguished by mode shape 
and natural frequency. For example, a case of flutter 
may involve the interaction of three modes: 
wing bending, one of wing torsion, and one of control 
surface rotation. It is further assumed that the action 
in any one mode is that of a simple mechanical system 
containing mass or inertia (the wing or control surface), 
linear springing (the elastic structure), and damping 
proportional to the amplitude of motion (the internal 
damping of the structure). The particular components 
making up the control surface system are the surface 


one of 
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itself and the structure connecting it to the command 
station. With the introduction of a servomechanism 
between the command station and the control surface, 
the mass-spring-damper idealization becomes untenable. 
The reaction of the servomechanism to excitation at 
its output end is a complex function of the frequency of 
excitation, a function that is best determined experi- 
mentally. In the present analysis, then, the idealiza- 
tion of the system governing control surface action is 
dropped. The system is allowed to have an arbitrary 
impedance (the vector division of sinusoidal exciting 
moment at the control surface by angular response of 
the surface). The conception of wing action in bending 
and torsion remains as in classical theory. 

Under these conditions, the equations of wing-con- 
trol surface motion are written. The method of solu- 
tion consists in choosing specific values of flutter fre- 
quency and flight velocity and determining the unique 
impedance that will allow the neutrally stable flutter 
condition to exist. In the solution, the flutter fre- 
quency and frequency of servo-oscillation are necessarily 
equal. The field of frequency and velocity of practical 
interest in a given case is covered by systematic calcu- 
lations of this sort. The values of the impedance re- 
quired for flutter are plotted as a set of loci, each locus 
representing possible flutter at a given frequency carried 
from zero velocity to the top speed of the aircraft. The 
loci are termed flutter boundary curves. Comparison 
with the experimental frequency variation of impedance 
of the servo and control structure allows the flutter char- 
acteristics of the aircraft to be deduced. 

The impedance of the classical control surface system 
is a theoretically defined quantity. Correlation with 
standard flutter theory is shown by consideration of 
this special case. The merits of the present method as 
applied to conventional (nonservo) aircraft are also dis- 
cussed. 

The technique that has been described is developed in 
detail below. Specifically, the work is divided into the 
following sections: (I) The Equations of Motion. (II) 
Specification of Flutter Boundary Curves: Case (a) 
Torsion-Bending-Control Surface Rotation; Case (b) 
Bending-Control Surface Rotation; Case (c) Torsion- 
Control Surface Rotation. (III) Sample Calculation. 


NOMENCLATURE 
Geometric Quantities 
a = distance from mid-chord to elastic axis in per cent of 
semichord, positive aft 
b = semichord 
c = distance from mid-chord to control surface leading 


edge in per cent of semichord, positive aft 
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Fic. 1. The wing-aileron system. 
= distance from mid-chord to control surface hinge line 
in per cent of semichord, positive aft 
h = vertical bending deflection, positive down 
a angular deflection in torsion of airfoil section, positive 
when trailing edge is depressed 
38 = angular deflection of control surface about the hinge 


line, positive when trailing edge is depressed 


Aerodynamic Quantities 


V = forward velocity of the aircraft 
p = density of the air 
Va = aerodynamic moment of the airfoil about the elastic 


axis per unit length of span at any instant in a steady 
oscillatory motion 

V3 = aerodynamic moment of the control surface about the 
hinge line per unit length of span at any instant in a 
steady oscillatory motion 

P = aerodynamic lift force per unit length of span at any 


instant in a steady oscillatory motion 


R,, = that part of the aerodynamic force or moment co- 
efficient produced in the rth degree of freedom by a 
unit displacement in steady oscillating motion in the 
sth degree of freedom that is in phase with the mo- 
tion in the sth degree of freedom 

I,, = the counterpart of R,, that is in phase with the velocity 


of motion in the sth degree of freedom 


Dynamic Quantities 


Ca = torsional stiffness of airfoil section per unit length of 
span 

C, = airfoil vertical bending stiffness per unit length of span 

Za = airfoil torsional damping coefficient 

g, = airfoil vertical bending damping coefficient 

I~ = mass moment of inertia of airfoil about elastic axis per 
unit length of span 

Jz = mass moment of inertia of control surface about the 
hinge line per unit length of span 

MJ’ = mass of airfoil per unit length of span 

Sa = mass static moment of the airfoil about the elastic 
axis per unit length of span 

Sg = mass static moment of the control surface about the 
hinge line per unit length of span 

We = natural frequency of a mode of torsional vibration of 
airfoil section 

w, = natural frequency of a mode of bending vibraticn of air- 
foil section 

w = frequency of flutter oscillation 

k = flutter parameter, wh/V 

t = time 

v = dv/dt, where v is a time variable 


(I) EgQuations oF MOTION 


Fig. 1 schematically portrays a wing-aileron system 


in displaced position. A unit span is considered and 
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the equations of motion for wing bending, torsion, and 
aileron rotation are written. The equations set down 
are not for an arbitrary motion but, as is customary, 
are for the specific condition of simple harmonic motion 
of all the variables. The condition is one of neutral 
stability, separating a region of stability from a region 
of self-excited instability or flutter. Consideration of 
the neutrally stable state allows the flutter boundary to 
be specified in terms of the parameters of the problem 


Summing the forces on the section: 


hM +h, + ig,) = P — &S. — BS l 


Complex notation is used to distinguish vector com 
ponents in phase and out of phase with the bending 
Idealized structural damping is indicated, 
such damping being in phase with the velocity of bend 


motion. 
ing and proportional to the amplitude. The damping 
constant is put in the form C,g, as a matter of con 
venience. 

Taking moments about the elastic center of the sec 
tion: 


al, (1 + 4.) = M. 


sa 


- i. 
Bl I, + ble — a) Sp] 2? 


ala + 


Lastly, taking moments on the aileron about its hinge 
line: 


(8B) = Msp — hSs — &[Ig + b(e — a)Sg] (3 


where the left-hand side of the equation is an unknown 
function of @. 
To formalize the condition of harmonic motion, let 


wl 7 
a= ave’ ] 
8 = Boe tA) | 
—— vy 
k= he" * af 
— Uv 


Substituting Eqs. (4) in Eqs. (1) and (2) and recog- 
oi { , / 7] 
nizing that w. = V C,/I, and w, = VC,/M, there are 


obtained 


yy 2 ; P 
- aSo —_ BSg a AM l— (“ (1 + 1g) — : (3 
~ Ww @” 
w.\* ; 
— ale} 1 — ( (1 + iga) | — 
WwW 
Me . 
8 Ip + db(e — a) Sg —hS, = - (6 
w? 


With the substitution of Eqs. (4) in the right-hand 
side of Eq. (3), the sinusoidal moment forcing the 
motion of the aileron system becomes 
Moment = Mg + aw*[Ig + b(e — a)Sg] + hw?Sg (7 
fg is harmonic and of the same frequency as 8.! 

The vector impedance of the aileron system may be 
defined in terms of its response to a sinusoidal forcing 


moment MN. Thus, 


Impedance 8 = IN/B s 


The impedance depends on frequency and, for the non 
linear servomechanism, on the magnitudes of IW and 
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3. The latter dependence is discussed at a later point Impedance 8 = 4; + 7 de (9) 
in the paper. Impedance £6 is expressed in real and Then 
imaginary parts as follows: M = B(I, + 7 Je) (10) 
Substituting in Eq. (10) the moment given by Eq. (7), the equation for aileron motion becomes 
8 , M 
—a[l, co ble — @& Se] + a (9; + 1 So) — hSg = = ( 11) 
w” w~ 
The aerodynamic moments and force are* 
— h Par ; 
Ma = —pbtw? | (Ray + tan) , + (Rag + t1a3)B + (Raa + tlaa) a 
ia , h ; ' 
Me = —apb ow” (Ror + tI) b + (Rog + lng) 8 + (Roya + 1 ne) , (12) 
. — iin eh ; | 
P = —7pb®w?| (Ra + tn) 5+ (Reg + Iea)B + (Rea + tex) 


where R,», Zan, etc., are functions of airfoil geometry and the parameter V/w. 


are available in both equation and tabular form.’ * 


Ren, Tan, etc., or their equivalents 


Substituting Eqs. (12) in Eqs. (5), (6), and (11), the flutter equations reduce to the following set: 


-. 
. apb* 
( Ig + b(e = 
a rob! 


ra 


a)S J So 
z ‘| + [Re + ilvel) + 3 ( Sak. 


w*arpb4 


— ra + (Re + il) + 


; a 
a i [Ree i ila} ) + 6 ¢ 
/ mpl 


mpb® 


h M w,\? ; 
- ~li— ) (1 + tg,) | + [Ro 
b pb w 


* The aerodynamic notation is that of reference 2. 


[ e (“*) ‘a+ ig) + i. + ile) + 2 (- | 


As derived, the equations relate to a wing section of — vanish. 
unit span. Following methods of conventional two- 
dimensional flutter analysis,‘ a sample wing strip is Q 
considered, and its action is assumed to indicate the C; 
action of the wing. The quantities we and w, take on Ce 
the meaning of the natural frequencies of the pair of C 
wing vibration modes that are being investigated for ee 
possible unstable combination with the aileron system. C. 
Extension to the more rigorous three-dimensional Ce 
treatment can be accomplished by integration of the C; 
unit strip effects along the span. Cs 

Eqs. (13), (14), and (15) are now solved for the values c 
of 9; and J. required for flutter. Cu 

Cu 
(Il) SPECIFICATION OF FLUTTER BOUNDARY CURVES Crs 
Cis 
Case\a) Torsion-Bending-Control Surface Rotation rao 

To satisfy Eqs. (13), (14), and (15), it is necessary C15 

that the determinant of the coefficients of the variables Ci 


Tp + b(e — a) Sg 


| + [Rig + ils) + 


mpb' 

h ka 
- -+ [Ran + tla) } = 0 (13) 

b mpb* 

+ [Rog + 7 hel) + 

h Sg . ; 
aaa : + [Rop + iD pp | = 0 (14) 

b mpb* 


(15) 


+ ifal) = (0) 


The following notation is introduced: 
l w- 

—(Iq/mpb*) + Roa 

Tea 

Tq Wa" /mpb' ° 
C280 

—{[Jp + b(e — a)Sg|/mpb'} + Ras 
lap 

—(S,/mpb*) + Ray, 

Lan 

— {t[Ig + b(e — a)Ss|/mpb*} + Roa 
La 

Rog 

Tyg 

— (Sg/mpb*) + Ro 

Ton 

—(S,/mpb*) + Rea 

Tce 
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Cz = —(Sp/rpb*) + Re Coo = Len 
Cis — Lg Cu = Me”, apb* 
Cio —_ —(M mpb’) + Re Coo = Cargn 4 


The equation in the determinant of the coefficients is then 


[C; + iC] 


QS; ) ( O45. 
C, 4 Lelc. 4 ™& 
( hid mpb* 7 =) 


1C\s| 


[(Ci + 2C:) + (Ce + 2C))) 
[Cy + 1C 0] 


[Cis i tC] [Cz 


It can be solved directly for 9; and Jz in terms of the constants C, and ©. 


imaginary parts, the following solution is obtained : 


J; = ( mpb A 


[C; 2 1C3] 





[Cis + tCy4] = 0 


(Cig + QC2) + 1 (Cao + QC 22) } 16) 


Expanding and separating real and 

















2) {[(4142 + BiB2)/(Ay + By)] — Cu} lv 


So = (mpb*/Q)}[(A1B., — A.B,)/(Ay? + By)] — Cp} IS) 


where 
A; = [CiCiy — C2Cxo — 
By = [CiCn + CC — GC — 


A, = [Ci(O3Gi7 — 


CsCig + Cris + CsCiz) + Cis(—CsCre + CoCr) + Cie(CsCu +. CeCizs)) + 2[C3(CisCr — 


CiC15 i CsCis] T 2(Ci C2 ae CoCr r C3Cig kone CxC29| + Q?[C3Co — CyCr2] 19 
CsCi5| + Q[ Ci Cx i CoCo + C3Coo + CxCi9] + 2?[C3Co2 oi C;C2| (20) 
CiyCis) —_ Co(CisCis + CuyCi7) + Co(C5Cig -— CoCr = CrCi7 T CsCis) ++ Cro( — CsCr nice 


CisCis) i 


Cy( Cia Cis + C14Ci7) + Cu (C5Co — CsCro) — Co2(C5Ci0 + C5Co) | Zz 


By. = [Ci(CsCis + CurCiz) + Co(CisCir — CuCis) 4 
C5Coy — 
Ca( C13 Ci; 


Case (b) Bending-Control Surface Rotation 


The torsional mode of vibration is eliminated by 
dropping the first row and column of the determinant 
in Eq. (16). 
bending, control surface flutter is given by Eqs. (17) 
and (18), with the following simplified values of the 
quantities A,, B,, A, and By: 


The resulting impedance required for 


A, = Co + 2Cx (23) 
By = Ca + QCxe (24) 
° Ay = C3Civ — CuCis (25) 
By = Ci3Cig + CuCi (26) 


Case (c) Torsion-Control Surface Rotation 


Tie last row and columa of the determinant in Eq. 
(16) is dropped to eliminate the bending mode. The 
values of A;, B,, As, and By pertaining to torsion- 
control surface flutter become 


A, = Ci + OC; 


B, = Cr + QC, (28) 
A, = C5Cy — CoCo (29) 
By = CsCww + Cel, (30) 


Co(CsCoo +- CeCi9 <a CrCig seas CsCi7) + Cio(CsCig piest 
CrCiz + CsCis) — Cis(CsCra + CeCis) + Cie(—CsCis + CoCis)] + 2(C3(CisCig + CiusCiz) + 
CiasCis) + Co1(C5Ci0 + C5Co) + Cx(C5Co ~~ CsCw) | (22 


Eqs. (17) and (18) give the value of the impedance re-_—f 
quired for flutter in terms of the basic parameters w 
V/wbd. 
craft-servo combination is explored by calculating loci 
of flutter at constant frequency, each locus being car 
ried from zero value of |’/wb to that value correspond 
The frequency 
range is taken from zero to a value high enough to pre 


and The entire field of operation of the air 


ing to the top speed of the aircraft. 


clude the coincidence of servo impedance and the im 
pedance required for flutter. 

A modification is made in the impedance before plot 
ting; the effect of the inertia of the control surface is 
removed. The impedance finally handled is that of a 
physical system identical to the true system but with 
an inertialess control surface. The inertia of the sur- 
face contributes an amount —J,/Q to J; and nothing 
to 2. The real part of the inertialess impedance is_ | 
then 9, + (/s/Q2). This procedure reduces the im 
pedance of the conventional mass, spring, damper con 
trol system to a quantity invariant with frequency t 
that is, a point on the flutter boundary chart. Thus, 
consider the equation of motion of the aforementioned 
svstem: ‘ 
8 

+ CyB = mM 3 


Ww 


168 + Cgge 
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FLUTTER OF 


where 


Cs = stiffness of the control surface in rotation 

gg = Structural damping coefficient for control 
surface motion 

sinusoidal forcing moment 


nu = 
The solution of Eq. (31) vields 


M/s = [(Ceg << (Ip Q2) | +7 (Cegp) (32) 


then, 


gs 
(33) 


= Cs — (lp) 
$2 = Cage f 


or, modifying 9,, the impedance components of the 
inertialess system become 


I, + (Ip 0) = Cal : 
. (34) 
I, = Cage f 
Eqs. (34) describe a point on the impedance plane. 
This condition reduces the solution of the conventional 
flutter problem to a simple case of impedance matching. 
The present technique applied to conventional flutter 
has the attribute of allowing a small band of critical 
frequencies to be examined for flutter with the param- 
eter I’/wb a secondary variable. It also makes im- 
mediately apparent, in a given example, the effects of 
varying control system spring and damping constants. 
This subject will be treated more fully in the section 
that follows. 

One last specification is made with regard to plotting 
It is most convenient to plot the 
that is, the logarithm of 
the magnitude versus phase angle. The two coordi- 
nates for the flutter boundary curve chart are then 


the impedance. 
logarithm of the impedance 


InV [$1 + (Jp/2)]? + 92? 
ind 
tan—'{92/[d1 + (Js/2)]} 


A sample problem is now used to expand and eluci- 

date the foregoing theory. 
(III) SAMPLE PROBLEM 

Consider an aircraft with the following character- 
istics of a chosen wing section: chord, 4 ft.; distance 
irom mid-chord to elastic axis, 0.67 ft. forward; dis- 
tance from mid-chord to leading edge of aileron, 1.072 
it. aft; distance from mid-chord to hinge line of aileron, 
1.200 ft. aft; 


|! = 1,447 slugs per ft. 

Sa = 0.941 slug ft. per ft. 
0.0645 slug ft. per ft. 
J, = 2.25 slug ft.? per ft. 
0.0394 slug ft.? per ft. 
155,000 ft.lbs. per rad. 


77,900 Ibs. per ft. . 


S = 


SERVO-CONTROLLED 





AIRCRAFT 





Za = 0 
a = 0 


It is desired to inspect the aircraft for flutter at sea 
level. Its top speed is 700 ft. per sec. at this altitude. 
Control system specification is given at a later point. 
The flutter boundary chart is prepared according to the 
following procedure: 

(1) The values of w for which loci will be cal- 
culated are chosen. They are 5, 10, 15, 20, 25, 30, 35, 
35, 35.5, 36.0, 40, 45, and 50 cycles per sec. Close 
spacing of the frequencies is taken in the region of the 
natural bending frequency. The range is terminated 
at 50 cycles. It becomes obvious, on plotting, that 
no flutter of practical interest occurs above this fre- 
quency with the particular servo and conventional con- 
trol systems involved. 

(2) 
frequency locus is determined by fixing a speed beyond 
which analysis is of no value. This speed is 700 ft. per 
However, the loci are car- 


The maximum value of V/wb for each constant 


sec. for the present case. 
ried to much higher speeds (uniformly to a value of 
V/wb of 8) in order to illustrate their mathematical 
behavior more clearly. The velocities involved in this 
extension invalidate the aerodynamic assumptions that 
These later portions of the loci, then, 
In other examples, 


have been made. 
have no physical significance. 
though, equivalent portions of the loci may occur in a 
valid speed range. The parts of the flutter boundary 
diagram of physical significance in this specific problem 
are differentiated on the diagram from the parts of 
mathematical interest only. 

(3) From the tabular values of the aerodynamic 
coefficients of reference 2 and the data of steps 1 and 2, 
a sufficient number of values of the constants C; through 
Cx» and A,, B,, Az and Bs are calculated to define the 
flutter boundary curves. No compressibility correction 
is used. 

(4) 
quency of 


The curves themselves, loci at constant fre- 


Inv [9 + ([g/2)]? + Io? 


versus 
tan l 1 Io [Si a (Ig Q)]} 
are plotted. 


Inv [5; + ([e/2)]? + 42? 
is represented by 
V Idi + (Ip/2))? + 52° 


plotted to a logarithmic scale. 

The flutter boundary chart is shown in Figs. 2a and 
2b. Fig. 2a deals with frequencies from 0 through 
30 cycles per sec.; Fig. 2b, with all higher frequencies. 
The loci are carried to 700 ft. per sec. as solid lines. 
Their mathematical extension is shown dotted. The 
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locus at 5 cycles per sec. can be calculated to only 500 
ft. per sec. because of the limit of V/wb equal to 8 in the 
tabular listings of reference 2. It is extrapolated to 
700 ft. per sec. The equations of reference 1 can be 
used to extend the aerodynamic tables if this is neces- 
sary. 

As an aid in the interpretation of the flutter chart, 
the values of impedance required for flutter at zero ve- 
locity are plotted as a function of frequency in Fig. 3. 
The figure traces the starting points of the flutter 
boundary curves. The phase angle of the impedance 
for the zero velocity case (physically trivial) is either 
0° or 180°, where the impedance is finite. Thus, re- 
ferring to Figs. 2a and 2b, the flutter boundary curves 
begin on the horizontal lines 0° and —180°. The zero 
frequency curve is a point at zero, infinitely far to the 
left on the 0° line. As the frequency rises, the starting 
point moves to the right until, at a little above 35 cycles, 
it is infinitely far to the right. In this remote region, 
the point shifts to the — 180° line and moves to the left 
with further increase in frequency. In the small range 
35+ to 37+ cycles (approximately the natural bending 
frequency), the starting point sweeps the —180° line 
from plus infinity to minus infinity. This cycle then 
repeats itself, centering this time about a frequency 
The 


starting point of the infinite frequency locus is at plus 


somewhat above the natural frequency in torsion. 


infinity on the 0° line. 

Fig. 2 in its entirety is complex. Consider its char- 
acteristics, first, from the point of view of conven- 
tional flutter. Allow the control system to take on 
definite values of Cg and gg. The impedance of this 
system is then a point on Fig. 2... The point, ordinarily, 
is close to, and above, the 0° line. The phase angle of 
the impedance, tan~'! gg, is usually less than 12°. For 
the sake of argument let gg = 0. The impedance is 
then a single point on the zero axis at a magnitude of 
Cs. For flutter to occur at some frequency, the flutter 
boundary curve at the frequency must coincide with 
Cg somewhere along its length. the 


value of Cg, then, flutter can only occur when a locus 


Regardless of 


crosses the zero axis at finite velocity or, in a limiting 
case, is tangent to the axis. There is no flutter, there- 
lore, at frequencies up to approximately 31 cycles. 
There can be flutter at frequencies above 31 cycles at 
10-cycles 
within the speed range of the aircraft. The frequency 
The mode involved 


some velocity but only between 35- and 
range of flutter is now established. 
is primarily bending, aileron. Each value of Cg has 
associated with it a definite flutter frequency and ve- 
The velocity is available through the known 
spacing of the quantity V/wb along each locus. Values 
of Cz from approximately 350 to 3,400 will cause flutter 
The complete 


locity. 


within the speed range of the aircraft. 
flutter specification of the ship can be deduced from 
Fig. 2 and drawn. Figs. 4 and 5 give the variation of 
flutter frequency and velocity with the natural fre- 


quency of the control system. Final results may be 
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stated as follows: the range of natural frequency of 
control system that will produce flutter within the speed 
of 700 ft. per sec. is 15 to 47 cycles per sec. The asso- 
ciated flutter frequency range is 34.5 to 39.5 cycles per 
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sec. The curves of Figs. 4 and 5 can be readily ex- 


panded into families for various amounts of damping. 

The existence of a higher frequency flutter mode is 
shown by the action of higher frequency loci that cross 
the range of Cg under consideration. The crossings 
are of no practical significance, since they take place 
at extremely high speed, but they illustrate the means 
by which two modes of flutter may be associated with 
a single value of Cg. 

To summarize the case of conventional flutter, it 1s 
seen that the technique of impedance plotting presents 
a rather powerful method of analysis. Where servo 
units are not involved, the loci may be computed only 
in the region of zero phase angle. Besides the evalua- 
tion of standard control systems, the flutter boundary 
plot can be used to investigate nonlinear spring and 
damping effects. 

The criterion for flutter in the case of servo-control 
is that the point representing servo impedance at a 
given frequency coincide with the flutter boundary 
curve of that frequency somewhere along the length of 
the The which detailed flutter 
boundary curves are required depend on the servo im- 
Critical frequency regions must be 
The flutter boundary curves that 


curve. regions it 
pedance path. 
closely inspected. 
begin on the 180° axis, for example, although covering a 
frequency range of only 2 odd cycles, sweep a large por 
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tion of the impedance plane. This complication is less 
severe in the conventional control case, for the curves 
become regular at 0° phase. 

Fig. 6 gives a response of a particular servomecha- 
nism with a quantity /s/Q added to its real part. With 
this characteristic, flutter will occur at approximate 
values of frequency and velocity of 36 cycles per sec 
and 700 ft. per sec., respectively. 

The servomechanism is a 
response will be different for different amplitudes of 
Of primary interest is the response at 
A response at some higher ampli- 
tude may indicate flutter. The flutter will 
however, only if an aileron disturbance causes the sys- 
In addition to 


nonlinear device. Its 
aileron motion. 


small amplitudes. 
occur, 


tem to pass this threshold amplitude. 
variation with amplitude, the servo response is usually 
a function of a number of electronic parameters. Ina 
given case, it is necessary to decide on the responses 
required from the conditions of the flutter investigation 
and from the particular servomechanism behavior. 
An experimental means must then be devised to meas 
ure truthfully these responses. The problems involved 
are jointly electronic and dynamic and are beyond the 
scope of this discussion. 

In summary, there has been presented a method of 
flutter analysis based on the concept of impedance 
mapping. It is hoped that broad usage will bring 
into proper perspective its relative merit with respect 


to other methods in the field. 


REFERENCES 


! Theodorsen, (Theodore, General The ory of Ae rody? 
Instability and the Mechanism of Flutter, N.A.C.A. T.R. No. 496, 
1935 

Jenscoter, L. K., and Libby, P. A., Tables for the Calculation 


of Aerodynamic Coefficients of an Oscillating Airfoil-Flap Systen 
in an Incompressible Fluid, Navaer SM-30, June, 1947. 
> Theodorsen, Theodore, and Garrick, 1. E., MWechan } 
Flutter, A Theoretical and Experimental Investigation of the Flutter 
Problem, N.A.C.A. T.R. No. 685, 1940 
‘Myklestad, N. O., Vibration Analysis; 
Company, Inc., New York, 1944. 


McGraw-Hill Book 





Notices should be sent directly to: 





Changes of Address 


Since the Post Office Department does not as a rule forward magazines to forwarding addresses 
it is important that the Institute be notified of changes in address 30 days in advance of publishing 


date to ensure receipt of every issue of the JOURNAL and REVIEW. 


Institute of the Aeronautical Sciences, Inc. 
2 East 64th Street, New York 21, New York 

















ientaaeded 


1 is less 


curves 


necha- 
With 
ximate 


er sec. 


e. Its 
des of 
nse at 
ampli- 
occur, 
1e sys- 
10n to 
tsually 

Ina 
ponses 
gation 
lavior. 
meas- 
volved 
nd the 


10d of 
dance 


bring 


espect 


ts) $96, 


1 Book 








6 eee 








Theories of Plastic Buckling 


S. B. BATDORF* 


National Advisory Committee for Aeronautics 


SUMMARY 


Che theory for the plastic buckling of columns, which appears 
finally to have achieved a satisfactory form, rests upon the well- 
established uniaxial stress-strain relation. The development of 
a correspondingly satisfactory theory for the plastic buckling of 
plates has been hampered by the nonexistence of an established 
polyaxial stress-strain relation in the plastic range. 

Present theories for the polyaxial stress-strain relation beyond 
the elastic range can be divided into two types, often called flow 
and deformation theories. Theories of plastic buckling based on 
deformation theories are in better agreement with experiment 
than those based on flow theories. On the other hand, tests in 
which a material is compressed into the plastic range and then 
subjected to shear at constant compressive stress are in better 
agreement with flow than with deformation theories. Legitimate 
doubt therefore has existed as to the validity of any theory for 
the plastic buckling of plates. 

As a result of studying these apparent contradictions, a new 
theory of plasticity has been developed which is of neither the 
flow nor the deformation type. It is based upon the concept of 
slip, and its formulation was guided more by physical, and less 
by mathematical, considerations than previous theories. Experi- 
mental evidence of limited scope but of crucial character is in 
better agreement with the new theory than with either flow or 
deformation theories. The new theory accounts for the apparent 
contradictions previously alluded to and justifies the use of de- 
formation theory in the analysis of the plastic buckling of plates. 


INTRODUCTION 


Mee THEORY FoR the plastic buckling of columns, 
which appears finally to have achieved a satisfac- 
tory form!~* rests upon the well-established uniaxial 
stress-strain relation. The recently established result 
that even a perfect column can buckle according to a 
tangent modulus law is rational and can be readily 
accepted on the basis of the following intuitive argu- 
inent. The buckling stress of a column depends upon 
its stiffness against bending. In the elastic range, the 
stiffness is the product of E, the elastic modulus, and J, 
the moment of inertia of the cross section. Beyond the 
clastic range, the stiffness depends upon the stress and 
also on the rate of increase or decrease of compressive 
stress accompanying the bending. The effective modu- 
lus at a given stress may be anything between the elastic 
modulus (which applies if the relative rate of unloading 
and bending is such that all portions of the cross section 
are unloading) and the tangent modulus Ey (which 
applies if the relative rate of loading and bending is 
such that all portions of the cross section are loading). 
Since a structure always buckles in the mode associated 
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with the lowest buckling stress and since the minimum 
buckling stress of the column is associated with the 
minimum bending stiffness, the effective modulus in 
plastic column buckling is the tangent modulus. 

The development of a correspondingly satisfactory 
theory for the plastic buckling of plates has been ham- 
pered by the nonexistence of an established polyaxial 
stress-strain relation in the plastic range. It is not 
that no theory exists but rather that theories exist in 
great abundance. 


PRESENT THEORIES OF PLASTICITY 


Mathematical theories of plasticity take as their 
point of departure the stress-strain law for simple ten- 
sion or compression. This is illustrated in Fig. 1. -In 
the elastic range, the relation between the stress o and 
the strain ¢€ is essentially linear and is expressed by the 
equation 

o=Ee (1) 
Beyond the elastic range, two methods of representation 
are used. In one, the finite relationship is retained be- 
tween stress and strain, but the elastic modulus F is 
replaced by the so-called secant modulus Es, which 
varies with the stress, so that Eq. (1) becomes 


Eg € (2) 


The other mode of representation is an incremental one 
employing the tangent modulus /;7, which also varies 
with the stress. The mathematical expression of this 


law is 


do = Ey de (:3) 





Stress 











Strain 


Fic. 1. Uniaxial stress-strain relation. 
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Eqs. (2) and (3) apply as long as the loading continues 
to increase. Unloading, however, occurs elastically 
that is, Eq. (3) applies with /, replaced by the elastic 
modulus F£, so that the relation between stress and strain 
during unloading is a straight line parallel to the elastic 
portion of the stress-strain curve. 

The various theories of plasticity seek to generalize 
these results so as to formulate correctly stress-strain 
relations under complex loading conditions in the region 
beyond the elastic range. Those that have been pro- 
posed up until the present time can be classified into 
two types: One generalizes Eq. (2) and therefore in 
volves relations of the finite type, while the other gen- 
eralizes Eq. (3) and therefore involves incremental 
relationships. The theories involving finite laws are 
sometimes called deformation theories; those involving 
infinitesimal laws, flow theories. Both types of theory 
assume that the plastic law proposed applies during 
loading, while unloading occurs elastically. Whether 
the change in the stress state constitutes loading or un- 
loading is supposed to depend upon whether the change 
increases or decreases a rotationally invariant function 
of ‘the stresses. This function, which is sometimes 
termed the ‘‘stress intensity,’’ is usually assumed to be 
proportional either to the maximum shear stress or to 
the so-called octahedral shear stress, which is related to 
the energy of distortion. In order to make a plasticity 
theory definite enough to be mathematically tractable, 
additional assumptions must be introduced, and these 
ordinarily include the assumption that the principal 
axes of stress coincide with the principal axes of the 
plastic part of the strain or its increment. (For a gen- 
eral survey of flow and deformation theories see refer 


ence 4.) 
THEORIES FOR PLASTIC BUCKLING OF PLATES 


It is only comparatively recently that attempts have 
been made to solve the problem of the plastic buckling 
of plates rigorously on the basis of one or another of the 
various mathematical theories of plasticity. A solution 
based on a flow theory was obtained by Handelman and 
Prager,® and a solution based on a deformation theory 
was obtained by Stowell.® 7 For the sake of simplicity, 
in the present paper attention will be limited to these 
two theories, which are reasonably representative of the 
types of plasticity theory on which they are based. 

The discrepancy between the predictions of the two 
buckling theories is, in many cases, marked, and the 
available experimental data have been in decidedly 
better agreement with the buckling theory based on 
deformation theory.» For example, Fig. 2, taken 
from reference 9, compares both theories with experi 
mental data on the buckling of long simply supported 
plates. On the other hand, theoretical objections have 
been raised to deformation theories as a class, and ex- 
periments have been proposed to determine whether 
deformation theories are in better 


flow theories or 
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Fic. 2. Plastic buckling of long simply supported plates 


agreement with experiment.'® '' Some experiments for 
this purpose were carried out by Roger W. Peters and 
Norris F. Dow at the Langley Laboratory of the 
N.A.C.A. In the most directly applicable of these 
experiments, a cylinder was compressed into the plastic 
region and then, holding the compressive stress con 
stant, was twisted. According to the deformation 
theory used by Stowell, which is closely related to the 
Nadai theory,'* the initial effective shear modulus 
should have been (/s//)G, whereas according to flow 
The 


experiment gave a value close to G in agreement with 


theories it should have been the elastic value G. 


flow theories and in disagreement with the deformation 
theory. With suspicion thus raised against the founda 
tion upon which is based the theory of plastic buckling 
in best agreement with experiment, there appeared to 
be grave doubts as to whether a satisfactory theory for 
the plastic buckling of plates could be said to exist at 
all. 

In some cases the evidence was so conflicting as to 
appear paradoxical. the theoretical 
buckling stress of an extremely long simply supported 
hinged flange in the elastic range can be shown to be 


For example, 


given by the formula 


G(t b)? | 


g¢= 


where ¢ and + are the thickness and width of the flange, 
respectively. Since the modulus appearing in this 
equation is that for shear, which is to be expected since 
the flange buckles by twisting, the formula should evi 
dently be amended in the plastic range by using the 
effective shear modulus of the material when subject 
to the compressive stress obtaining at the moment of 
buckling. This effective shear modulus was found in 
buckling experiments to be (/s/ £)G, in agreement with 
the predictions of the deformation theory (more pre 
sisely, deformation theory gives (Es/E)G times the 
nearly unit factor 2(1 + w)/[3 + (2u — 1)(Es/E)], where 
nis the elastic value of Poisson’s ratio; for simplicity, 


ry 
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this factor will be ignored throughout the present 
paper). On the other hand, the effective shear modulus 
determined directly by compressing and then twisting 
the cylinder, as just mentioned, was found to be the 
elastic value G in agreement with flow theory. 


A New THEorY OF PLASTICITY 


In delving into the reasons for this paradox, a new 
theory of plasticity was developed which is of neither 
the flow nor the deformation type.'* This new theory 
assumes that all plastic deformation is due to slip, 
which occurs first along the plane of maximum shear 
stress in the direction of this stress at the moment the 
elastic limit for the material is exceeded. As the stress 
continues to increase, the slip along this plane also 
increases. In addition, slip commences along new 
planes slightly inclined to the plane of maximum shear 
as the shear stress on them also goes above the elastic 
limit. Moreover, as the shear stress on a given plane 
passes above the elastic limit, the number of directions 
in the plane in which the component of shear stress is 
large enough to cause slip increases. The theory as- 
sumes that slip occurs along every plane and in every 
direction in these planes in which the corresponding 
shear stress component is large enough to cause slip. 
The amount of slip in any direction along any plane is 
assumed to depend only on the prior history of the 
component in the direction of slip of the shear stress 
on this plane. Thus, strain hardening is anisotropic; it 
will take place only in those directions and planes along 
which slip has occurred. 

The anisotropic character of strain hardening in the 
new theory is nicely illustrated in the example of a 
specimen subjected to a loading history in which pure 
tension is initially in a given direction after which the 
direction of the pure tension is rotated, keeping the 
stress nagnitude constant. (Such an experiment was 
proposed for the purpose of assessing the relative merits 
of flow and deformation theories in reference 11.) Ac- 
cording to flow theories, there is no loading during the 
rotation of the axes because there is no change in any 
stress invariant, and there will therefore be no addi- 
tional plastic strain. According to deformation theor- 
ies, the situation is ambiguous,'* '° since a finite rota- 
tion accompanied by an infinitesimal increase in stress 
magnitude would give a finite plastic strain, whereas 
the same rotation accompanied by an infinitesimal 
decrease in stress magnitude would give a zero plastic 
strain, since unloading occurs elastically. What is im- 
portant for the present purpose, however, is that in the 
latter case both types of theory, and indeed any theory 
that defines loading in terms of the increase of a stress 
invariant, predict zero plastic strain during rotation. 
The new theory, on the other hand, predicts finite 
plastic strain, for during the rotation of axes new and 
unhardened planes will acquire shear stresses above the 
elastic limit and will therefore slip. 
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Fic. 3. Compressive stress-strain relations for cylinder subjected 
to compression and torsion. 


An experiment somewhat like the one just described 
has been performed at Langley. A cylinder was com- 
pressed beyond the elastic range and was then twisted 
while holding the compressive strain constant. The 
essential data are given in Fig. 3. Throughout the 
twisting process (AB in the figure), the compressive 
stresses decreased while the shear stress increased. 
The compressive plastic strain continually increased. 
Simple calculations show, however, that during the 
early part of the twist both the octahedral shear stress 
and the maximum shear stress were decreasing. Com- 
patibility with the known uniaxial behavior implies 
that, when the stress intensity decreases, plastic action 
stops and is not resumed until the stress intensity ex- 
ceeds the highest previous value. If, as in most plas- 
ticity theories, the stress intensity is assumed to be 
proportional to the octahedral shear stress, the value 
at the beginning of twisting was not exceeded until the 
shear stress exceeded 8.6 ksi., and, according to these 
theories, the compressive stress and strain during twist 
should therefore have been related by the line AC. 
Maximum shear stress theories would predict the same 
relationship up to a shear stress equal to about 5 ksi. 
No attempt has been made to calculate the subsequent 
behavior on the basis of flow and deformation theories 
because there seems to be no general agreement as to 
how to apply these theories beyond the first loading and 
unloading. The predictiors of the new theory as to 
the compressive strain associated with the applied 
stresses are in far better agreement with experiment 


(see Fig. 3). 


BUCKLING ON Basis OF NEW THEORY 


The new theory affords an explanation to the paradox 
discussed earlier. An ideal long hinged flange loaded 
in compression remains straight until the moment of 
buckling, when it commences to twist. The twist may 
take place at constant compressive stress or may be 
accompanied by an increase in compressive stress 
(shown as type A and type B changes, respectively, 
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in Fig. 4). According to the deformation theory, the 
effective shear modulus is in both cases (Es/E)G. Ac- 
cording to flow theory, the effective shear modulus is in 
both casesG. The new theory, on the other hand, gives 
an effective shear modulus that depends upon the type of 
change under consideration. According to the new the- 
ory, many of the planes that were slipping when the ap- 
plied stress was compressive commence to unload when 
shear is added in a type A change. Moreover, most of 
the planes that continue to slip are so oriented that they 
make only a small contribution to plastic shear deforma 

tion in the original system of coordinates. The effective 
shear modulus is therefore close to the elastic value G; 
just how close depends on the compressive stress prior 
to twisting. For a stress increment of type B ap- 
proaching the vertical, on the other hand, practically 
all of the planes continue loading and suffer the same 
plastic deformation as though the final state had been 
reached by keeping the stress ratio fixed at the final 
value throughout the loading history. For loading 
conditions in which stress ratios are kept fixed, deforma 

tion and flow theories predict the same results, and, for 
the case in point, both give (/s/E)G for the effective 
shear modulus. Therefore, in the new theory, a type 
A change results in an effective shear modulus approxi 

mately G, and a type B change is close to the vertical 
results in an effective shear modulus approximately 
equal to (Zs/E)G. Since the twist experiment first 
discussed was a type A change and buckling is, accord 

ing to the argument in the ‘Introduction’ to this 
paper, a type B change, the new theory accounts for 
the apparent contradiction and also justifies the use of 
deformation theory in the analysis of the plastic buck- 


ling of plates. 
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Note on the Relations 


Between Viscous and 


Structural Damping Coefficients 


WALTER W. SOROKA* 


University of California at Berkeley 


SUMMARY 


The use of the structural damping factor g in vibration prob 
lems implies that the damped free vibration frequency of the 
structure increases with damping. Relations are derived com- 
paring the natural frequencies and logarithmic decrements for 
viscous and structural damping. For low values of damping (up 
tog = 0.20), the two concepts give nearly the same results when 
cis taken as twice the critical viscous damping ratio. 


INTRODUCTION 


b iben MECHANISM whereby an oscillating system 
dissipates its energy of vibration is one of the most 
important properties of the system. Whether or not 
it is desirable, it is nonetheless present in all mechanical 
systems. It is certainly useful in those cases where it 
limits the resonant amplitudes of structures and drive 
systems; where it stabilizes instruments, servomecha- 
nisms, and automatic control systems; and where 
it reduces the danger of flutter. In such cases, con- 
trolled amounts of damping may be deliberately intro- 
duced when inherent damping is insufficient. 

In general, analytical discussions of damping in vi- 
brating systems are concerned with damping forces 
proportional to velocity—i.e., proportional to the prod- 
uct of frequency by amplitude. This has been par- 
ticularly true of problems connected with engine and 
instrument isolation, torsional vibration of crankshafts, 
ete. More recently, there has been introduced the con- 
cept of a damping force proportional to amplitude! ? 
but independent of frequency. Experience has indi- 
cated that, in structures such as aircraft wings and 
fuselages, the damping loss depends more closely on the 
latter relation. In both cases the damping force is in 
phase with the velocity of vibration. 

Little has been written regarding the implications of 
this concept of damping in structures and how it com- 
pares fundamentally with the concept of viscous 
damping. It is the purpose of this note to compare 
the newer treatment with the classical method and to 
bring out the main points of difference and similarity. 


FREQUENCIES OF THE TWo SYSTEMS 


The classical treatment for the free vibrations of a 
single degree of freedom damped system sets up the 
differential equation of motion in the form 
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mx + c&% + kv = 0 (1) 


where m is the mass of the system, ¢ is the viscous 
damping coefficient, and & is the spring stiffness. 

The treatment propounded in reference 1 sets up the 
equation of motion in the form 


mk + k(1 + ig)x = 0 (2) 


where g is the dimensionless structural damping co- 
efficient andi = V —1. Here a complex stiffness, 
k(1 + ig), serves to represent both components of the 
force proportional to displacement, x. One compo- 
nent, kx, isin phase with the amplitude; this is the usual 
spring force. The other component, ikgx, is in phase 
with the velocity (90° ahead of the displacement) by 
virtue of the unit vector 7 introduced into Eq. (2); this 
is the damping force. 

The solution of Eq. (1) is presented in detail by 
Den Hartog.* This solution is 


x = e~/?™(C, cos w,tf + Co sin w,f) (3) 


wherein w, is the natural frequency of the damped 


motion and is given by 
on = Vk, m) — (c?/4m?) (4) 


The damping coefficient may be put in dimensionless 
form as a ratio taken with respect to the critical damp- 


: ah / n : 
ing coefficient, c. = 2V mk. Thus, with 
A = ¢/C; (>) 
/ ' . 
wn, = VkR/m Vi — »? (6) 


The damped natural frequency decreases, of course 
as the damping coefficient increases. 
To solve Eq. (2), let 


v = Xe" (7) 
Substitution of Eq. (7) into Eq. (2) yields 
ms? + k(1 + ig) = 0 (8) 


In Eq. (8) either g is zero or s must be complex to 
satisfy this relation. Since g is a property of the sys- 
tem, the latter situation must prevail. Therefore, let 

s=a+t+ib (9) 
where a and 6 are to be determined so as to satisfy 
Eq. (8). Then, 


m(a? + 12ab — b*) + k(1 + 4g) = 0 
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which splits into the two equations: 


2 abm + kg = 0) 


: (10) 
m(a*? — b*) +k = of 
From Eqs. (10) we derive the fact that 
a? = (k/2m)(-1 + V1 + g?) (11) 
b? = (k/2m)(1 + V1 + g?) (12) 


In Eq. (11), the negative value of a is to be taken, 
since this value represents a decay of the disturbance 
whereas a positive value of a represents increasing am- 
plitude, which is a physical impossibility in a damped 
free oscillation. 

Eq. (12) is an equation for the natural frequency of 


the damped free motion. We may rewrite it: 


eo, = Vk 2m V1 + Vi + g” 


A fundamental difference between the two assump- 
tions of the character of the damping becomes evident 
from a comparison of Eqs. (6) and (13). Whereas the 
assumption of damping proportional to velocity implies 
a natural damped frequency that decreases with increase 
in damping, the assumption of damping proportional 
to amplitude implies a natural damped frequency that 
increases with increase in damping. The latter 
phenomenon is certainly not consistent with the usual 
conception of the effect of damping on a free vibration. 


(13) 


LOGARITHMIC DECREMENT 


The envelope of the free vibration amplitude is given 


wot / : ~ 
by e ““*, where w = Vk/m for the first case [see 


tq. (1)]. 
pression is é€ 
a free vibration is afforded by the ratio of amplitudes 


For the second case, the corresponding ex- 
— . A measure of the damping in 


at the beginning and end of a complete cycle. The 
time interval for a cycle is given by 27/w, sec. Thus, 
for the first case, the amplitude ratio, R, would be 
Awol 
R = . po hwo /wr p2ad/V/1—X 14) 
- seni € ( 
e 
In the second case, 
R=em™ e2ra/(l+ V1+e") 15 
The logarithmic decrement, 6, in each case is 
6 = 2r\/V 1 — DX? (16) 
6 = 2rg/(1 + V1 + g?) (17) 


In general, A is small in engineering structures, per- 
haps of the order of 0.02-0.04 or less, when no effort is 
made to introduce damping deliberately. This makes 
— )?* little different from unity. 


In flutter analysis the 


the quantity V 1 
Similarly, g is extremely small. 
value of g is rarely taken in excess of 0.05 and usually is 
this makes the quantity 


appreciably less. Again, 
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/ eine, ; é 4 
1+ V1 + g? little different from the value 2. Con- 
sequently, for small damping, Eqs. (16) and (17) may 
be rewritten, approximately, 


6 = 


27 X 
6 = xg 


Whatever assumptions are made for the damping 
effect in any given system, the rate of decay of a free 
vibration on test will be a physical property of the sys- 
tem irrespective of the assumptions. Consequently, 
the approximate relationship between A and g for small 
damping from Eqs. (16a) and (17a) is given by 


g = 2d (18) 


Table 1 uses the relationship (18) to provide corre- 
sponding values of \ and g for comparing the logarith- 
mic decrements and damped natural frequencies of a 
structure under the two assumptions of damping. 


TABLE | 
Comparison of Logarithmic Decrements and Natural Frequencies 


Viscous Damping Solid Friction Damping 


ny 6 Wn/ wo g(= 2A) 6 wn/w 
0.01 0.0628 1.00 0.02 0.0628 1.00 
0.02 0.126 1.00 0.04 0.126 1.00 
0.05 0.315 1.00 0.10 0.3138 1.00 
0.10 0.631 0.995 0.20 0.622 1.00 
0.20 1.28 0.980 0.40 | 21 1 ow 
0.50 3.63 0.866 1.00 2.60 1.10 
0.80 14.0 0.600 1.60 3.48 1.20 


HYSTERESIS Loops 


When a structure is subjected to steady-state vibra 
tion, the tip of the resu‘tant restoring force vector (com- 
posed of the spring force and damping force) traces 
out an ellipse when plotted against displacement. This 
is true for either assumption as to the form of damping. 
The equations for the ellipses are given in parametric 


form by 


F = xok[sin wt + (cw/k) cos wt] | 19 
. 4 \ : 
x = Xo sin at f 
and 
F = Xok(sin wt i £ Cos wt) ( (P0 
x = Xo sin wi f - 


where w is the steady-state forced frequency and x9 1s 
the amplitude of the motion. 


The areas of the ellipses will be given by 


A= $F dx = rCwx,’ 21) 


A = tkg Xo" (22) 

Eqs. (21) and (22) represent the amounts of energy 
dissipated per cycle under the two assumptions of 
damping. As the frequency w goes to zero, the ellipse 
given by Eq. (19) gets progressively thinner, reaching 


zero at the same instant that w reaches zero. In the 


(Continued on page 448) 
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ABSTRACT 

The problem of interference between wing and body which was 
investigated in the first report of this series (JOURNAL OF THE 
AERONAUTICAL SCIENCES, Vol. 15, No. 6, p. 317, June, 1948) is 
again examined, but now, in contrast to what was done there, 
the small perturbation condition is eliminated—i.e., it is no 
longer required that the perturbations produced by the body and 
the wing in the uniform flow field be small enough to permit the 
linearization of the differential equation of the fluid motion. It 
is still assumed, however, that a potential field exists—a condi- 
tion less restrictive than the previous one—and the determina 
tion of the flow is carried out under the assumption that the very 
tip of the nose is conical. 

Since the free-stream velocity is supersonic, the field in the 
region ahead of the wing does not depend upon the presence of 
this appendage, and thus this phase of the study warrants the use 
of the method already propounded by the author in a previous 
report, to which, however, substantial modification and simplifi 
cation are now applied. As in the first report of this series, the 
wing is assumed to extend beyond the region perturbed by the 
body, or at least it is assumed that the wing plan form is such that 
the tip phenomena due to the finite aspect ratio of the wing re 
main unaltered in presence of the body. 

The determination of the fluid motion is accomplished through 
use of the hodograph method in derivation of the field correspond 
ing to the conical nose and through application of the method of 
characteristics in the remaining regions. 


INTRODUCTION 


F  gumeresage THE EFFECT of interference between wing 
and body at supersonic speeds was treated with a 
satisfactory degree of precision in Part I of these reports, 
wherein the method of small perturbations was em- 
ployed, the restrictions of that approach may not be 
tolerable except in cases where the results are required 
with a relatively small amount of labor and where ac- 
curacy and generality can be sacrificed to expediency. 
A more exact solution to the problem is afforded by 
the method of characteristics, which is expounded in 
the present report so as to shape the details of this 
fundamental mathematical procedure to the specific 
three-dimensional problem requiring analysis. The 
report will thus have the nature of a stepwise procedural 
development of essentially simple relationships that 


are used repetitiously; consequently, the work will 
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appear lengthy, while in fact there is little of a complex 
character involved. Most attention will be devoted 
to the working out of a systematic and thorough cover 
ing of every region of flow by a network of characteristic 
surfaces, from which the local flow conditions may then 
be ascertained. 

Even this general and precise investigation of the 
problem of wing-body interference is, of course, based 
on several assumptions that may not agree rigidly with 
physical reality. An important restriction imposed 
on the present method is that the flow over the body is 
irrotational even in the presence of the shock wave 
emanating from the nose. It is believed that the dis 
crepancy between physical reality and the present 
theory will be sufficiently small in most instances that 
it will not warrant the development of the more intri 
cate solutions for which a potential does not exist. 

Further restrictions to the method as applied here are 
artificial in that they can be removed or altered so as to 
fit any particular example that does not conform to the 
configuration being treated. Of course, the body has 
been assumed axially symmetric and also unyawed, 
although the angle of attack and wing incidence are 
arbitrary. The body contour, if not made up of jointed 
straight-line segments, is approximated thusly, so that 
the body is composed of a conical tip and a series of 
frustums of cones for the nose, while farther back a 
cylindrical portion is finally reached in the neighborhood 
of the wing-body juncture. 

Because the flow is supersonic, the presence of the 
wing does not affect the fore part of the body, and, 
consequently, it is possible to begin the investigation 
of the entire flow field surrounding the body by means 
of standard procedures, which will be repeated here 
first of all for the sake of completeness. 


List OF SYMBOLS 


Vi = limiting velocity of the stream (velccity at 
tained upon adiabatic expansion of the flow 
to an absolute vacuum) 


Vo = velocity of the undisturbed stream 

C = velocity of sound in the undisturbed stream 
G = local velocity of sound in the disturbed stream 
vo = V/V 

co = O/V; 

¢ = Gy V, 








412 





Fic. 


1. 


JOURNAL 


OF THE 


Orientation of rectangular coordinate system. 





| 
| 
| 


_————— 





——— = 


‘46. 2. Coordinate and vectorial orientation in the spherical 


(x, y, 0 


coordinate system. 


the potential of the field around the conical 
nose 

the potential for the field that exists down 
stream of the initial characteristic surface 
from which perturbations are propagated 
into the fluid in consequence of deviations 
in the shape of the body from the conical 
form 

the velocity potential in the region down- 
stream of the initial characteristic surface 
(or shock wave) from which perturbations 
originate because of the presence of the wing 

Cartesian coordinates with x along the axis 
of the body and oriented (Fig. 1) in the 

the undisturbed velocity 


same sense as 


> 
Vo; the Z-axis is normal to the x-axis and 
in the plane of symmetry, oriented in the 
sense of positive lift; the Y-axis is normal 
to the plane of symmetry and directed ac- 
cording to a right-handed system; the 
origin is located in such a way that the lead- 
ing edge of the wing lies along the Y-axis 


cylindrical coordinates with x situated as for 


the Cartesian coordinates; the initial side 
of 6 is the direction of the positive Y-axis; 


and y is perpendicular to x in any meridian 
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plane and is positive in the direction of the 
terminal side of 6; the origin is at the tip of 
the conical nose (see Figs. 1 and 5) 
(r, ¢, 0) = spherical coordinates with origin at the apex 
the initial 
side of the angle ¢ is in the direction of the 
positive x-axis (see Fig. 2) 
velocity components along the Cartesian axes 
velocity components based on the cylindrical 
coordinate system 


of the conical tip of the nose; 
(Vz, Vy, Vz) = 
(Vz, Vy, Ve) = 


velocity components based on the spherical 
coordinate system (Fig. 2) 


(Ves Vo» Vo) = 


The same letters in lower case refer to the ratios of the corre 
sponding velocities to the limiting velocity V,, while the use of 
primes or various subscripts will always denote velocity compo- 
nents at the points indicated in the text. 


> 


By = angle of inclination of Vp with respect to the 
x-axis 

p = pressure 

p = density 

k = adiabatic exponent (taken as 1.400 for air) 


The meaning of the other symbols not listed here is indicated 


at appropriate points in the text. 


(A) DETERMINATION OF THE FIELD AROUND 


THE CONICAL NOSE 


(A-1) Potential of the Conical Field 


(A-1-1) General Equations. Let the potential of the 


conical field be expressed as 

®,V, = VirF(¢, 8) 
The 
velocity components 2,, v,, v, along the radius 7, in the 


in terms of the spherical coordinates (7, ¢, 6). 


direction normal to r and in the meridian plane, and in 
the direction normal to the meridian plane, respec 
tively, are given by: 


ss Oo; — uh l OP, OF 

or ~ +e Oe 
sds OhtiL:COOOF é 
v“oeo=— = 


rsing O86 sin ¢ O60 


The differential equation defining ®,; is, consequently 
i 5 : 


— VCOr tt, OF .. 
sin? ¢ { | t - 2— sin ¢ + 
= Oy" c” O86 Oy 
(1 ; r) oF ee (2 _ & + ") 4 
cy a c 


= Vor 9 
a sin 2y (1 + * | = 0) 
“ c~ 


while the local velocity of sound, c, is related to the 
velocity of sound in the undisturbed stream through 
the Bernoulli equation: 
' oa Bork : 2 ; 
cae) = & ~ 6 = U9”) 
(A-1-2) Determination of the Function F by the Hodo- 
graph Method. 


If, in every meridian plane @ = con- 
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INTERFERENCE 
stant, the hodograph of the motion is drawn in terms 
of the components v,, v, of the velocity in the direction 
of the cylindrical coordinates x and y, it is evident that 
the velocity of the fluid for all points of the polar radius 
m, which is inclined to the x-axis by the angle ¢ (Fig. 
3), is given by the tangent to the hodograph at the 
point M, and this tangent is perpendicular to the radius 
itself; also, the radius of curvature of the hodograph 
at M is:? 


R = v, + (0v,/0¢) (4) 


Because of Eq. (4) and taking into account Eqs. 
(2),"it is possible to write Eq. (3) in the form: 


UV," 2 v, oe rt 
(1-")r=—' at (Zr) — cot x 
c* sing c*? \ 06 


U9" U9" 1 oy 
(1 +— |-(: _ “)(», a ~*) 3°) 
E ( i =) - - sin ¢ 00 


This equation permits the determination of the hodo- 
graph in all the meridian planes if the values of the 
velocity components v,, v,, v are known for all the 
points of a given surface S, which is assumed to be an 
arbitrary cone with apex at the origin. 

Let the equation of this surface be: 





¢ = ¢(9) (5) 
and let the known values of vy and v, on S be denoted 
by v% = u(6); v1 = 2,(6); then 


dvp 4. 2% 
(F) = 20 *3. ¢ ) 
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Graphical extension of the hodograph to successive 
reference cones for flow over the conical nose. 


where ¢ is the derivative with respect to @ of the func- 
tion that is defined by Eq. (5). 
But from Eq. (2): 
Ov, 1 Ov, 


Oy sin ¢ O06 


, 1 ov 
¢|— “—cot¢ | 
sin ¢ O86 


— cot ¢ (%) 
and, hence, 


(“*) _ de 5 
d6 Ss 06 


Analogously, 


(). a ae 
£ = (R-y», 
i: = 6h 6S vy 


whence, 


dv,/d0 = (dv, 


and therefore, with aid of Eq. (3’) 


V_" Ss 2 Ye [1 — (v?/c?)] ., 
£ e+ 71 g? 


- = — ¢ _ 
c sing ¢ sin? o sin ¢ 


'd0); — (R —1,)¢ 


¢g dv, . 
.  ( +) + vu¢ — cos¢ co | 
y Ss 


O% R +(@) 
00 ‘ sin ¢ dé/s5 ~~ sin 


VV dv, V9" 
2[ (tp), + ne] -oeora(+2)- 
c* dé Ss c* 


ve?\ { 1 (dvs ¢ dv, r 
we ) $2, ( A sa" (“ + v,¢ — cos ¢ (u%) | ¢ 6) 
( c? \ + sin g \ dé sin? mI d6/ 5 ° sities 


The second member of Eq. (6) is now known for all 
values of 6 and for the value of ¢ which corresponds to 
the intersection m of the above-defined conical surface 
S with the meridian planes. Having located, for each 
meridian plane, the point M, which represents the 
velocity corresponding to the straight line m (Fig. 3), 
it is at once possible to locate the center of curvature, 
C, of the hodograph at M, since it lies on the parallel 


to m through M at a distance R, which may be calcu- 
lated from Eq. (6). 

Thus it is possible to draw, 
planes relative to the various planes @ = constant, the 
arcs of circles that are osculatory to these hodographs, 
between the straight lines ¢ and ¢ + Ag (Fig. 3). 
Thus, the values of the velocity components v, and 
v, are determined for all the straight lines m’, inclined 


for all the hodograph 
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Fic. 4. Arrangement of the shock surface at the nose 


at the angle g + Ag, while the values of vj are obtained 


from the equation : 


Ov, 
Va (Uy y=¢@ + Ag = (Vale ¢' | cot yay) {- 
Oy e~*¢ 


l (“*) SR | : ” 
D\ a v,) pa ¥ts) (¢) 
sin ¢ d0/ s aj e=¢ 


In this manner the values of v,, v7, and v,, relative 
to a new conical surface S’ defined by the equation 
¢ yg + Ag, are obtained. The application of the 
above-described procedure can thus be repeated, and 
it is apparent that one can continue the drawing of the 
hodograph relative to all the meridian planes by pro- 
ceeding step by step. If ¢ is extremely small for the 
surface SS (i.e., if S closely approximates a right circular 
conic surface) and under the assumption that Ag is 
taken to be different from one meridian plane to the 
next, the problem reduces simply to that of a right 
circular conic surface as far as S’ (or the successive 
S’’) is concerned; this case makes possible an immedi- 


ate simplification of Eq. (6). 


(A-2) Shock Wave on the Conical Nose 


(A-2-0) Alternative Points of View. The problem of 
determining the velocity field around the conical nose 
is now reduced to that of establishing the initial data 

that is, v,, v, and von a conical surface S. Two 
procedures are available for this purpose. The first 
consists of assuming the data appertaining to the un- 
that is, the ratio Vo/V; and the 
angle of attack 6,—and then in determining the shock 
waves that are compatible with such a stream and with 
a conical field of the sort dealt with in (A-1). The 
surface possessing the initial data is then the shock sur- 
face itself, while the form of the conical nose of the 
body is consequently determined as a result. The 
second procedure consists of assuming the conical shape 
of the body, a priori, and also the initial conditions of 
the stream over it; the surface possessing the initial 


disturbed stream 
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data is then the surface of the body itself, and the ar- 
rangement of the shock wave, as well as the condition 
of the undisturbed stream, is consequently determined. 
The first procedure just mentioned is developed below 
in Sections (A-2-1) and (A-2-2) and is valid for a general 
conical shape of the nose of the body with the sole limi 
tation that the surface be symmetrical with respect 
to the plane 6 = x/2. In contrast, in Section (A-2-3), 
the second procedure is dealt with, but there it will be 
necessary to restrict ourselves to the case in which the 
nose 1s part of a right circular cone with its axis coinci 
dent with the x-axis. 

(A-2-1) Possible Conical Shock Waves in Uniform 
Flows. — Let 


g = (8) Ss 


be the surface =, which is the shock wave front. Let 
mM, M2, m3 denote the straight lines that are, respec 
tively, the intersection of = with the general meridian 
plane 6 = constant, the perpendicular to m, in the 
meridian plane, and the perpendicular to the meridian 
plane (Fig. 4). Then the components of the undis 


turbed velocity v along ™, m2, m3 are: 


For m,: % (cos 8; cos ¢ + sin B; sin ¢ sin 6) 
For m2: wv (—cos B; sin g + sin B; cos ¢ sin 6) 


For m3: vo sin B; cos 8 

while the analogous components of velocity after the 
shock are: 

For m: F 


OF Oy 
(1/sin ¢)(OF/06) 


For mz: 
For m3: 


Imposing the condition, therefore, that the change in 
velocity across = be perpendicular to = itself and 
knowing that the direction cosines of the normal, 7, 
to = at a general arbitrary point are (with respect to 
Mm, Ms, and m3 which you pass through the same 


point) : 
cos(n, m,) = O 
cos(”, m2) = sin ¢g V sin? g+f? 
.  — = wa 
cos(n, m3) = —f/V sin’? ¢ + f? 
where f = df/d6, then 
v(cos 8; cos g + sin 6; sin g sin 8) — F 
0 
vo(—cos B; sin g + sin 6; sin 6 cos ¢) — (OF/d¢) 
sin ¢ 


v sin 6; cos 6 — (1/sin ¢)(O0F/08) H 
—f 
where H is the common value of the indicated 


ratios. Then on =—that is, for ¢ = f(@)—it follows 


that: 


v 
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—" F = v(cos 8; cos g + sin f; sin ¢ sin 6) It can be verified directly that Eqs. (9) are compatible 
ition dF/de = w(—cos B, sin g + (9) because the expression of F, OF /0@, and OF/O0¢g by 
a sin 6, sin 8cos g) — Hsing ‘their equivalents from Eqs. (9) confirms that dF/d@ = 
selieies vF/d00 = vw sin B, cos @ sin g + Hf sin ¢ (OF/00) + (OF/Od¢)(f). 
neral 
limi 
spect 
2.9). On the other hand, the equation for the shock must be satisfied through =. Let v,; and v,2 be the velocity com- 
ill be ponents (divided by V;) along the normal to = of the flow before and after the shock, respectively; and let 2, 
1 the be the component (also divided by V,) of the velocity which is tangent to the shock and which remains unaltered. 
INCI Then the equation of the shock has the form: 
UniUne = [(R — 1)/(R + 1) JQ — 2,/*) (10) 
form 
Now 
é 
i 
8) | ; - sin ¢ : 
Uni = Uo(—cos f; sin ¢ 4- sin #; sin 8 cos g) —j—————> — sin B, cos 8 => 
Let Vsin?g +f? V sin? gt+f? 
ad and 
dian 7 ; 
tie a( sin ¢ ) 1 OF ( f ) 
: V2 = >|  ] - = 7 aaa 
dian O¢\V sin? g+f sin ¢ 06 V sin? g +f? 
idis . 
while 
; 2 Jy [sin g (— cos 8; sin g + sin B; sin 8 cos ¢) — f sin 8; cos 6}*) 
Vv = Uo* — 7 > 
: sin’ g + f* f 
In addition, let 
a = sin ¢ (sin 8; sin 6 cos ¢ — cos B; sin g) — sin B, f cos 6 (11) 
the 
Then, substituting these expressions for v,, ¥,2, and v, into Eq. (10), there evolves: 
H = [2/(k — 1)] [avo/(sin? ¢ + f?)] — [(k — 1)/(k + 1)][C1 — m?)/av] (12) 
which, together with Eqs. (9), defines all the postshock velocities that are permissible as a result of shock action 
(upon the uniform stream under consideration) by agency of the conical waves [Eq. (8) ]. 
e im 
and =f 
1m (A-2-2) Complete Determination of the Conical Field normal to the contour surface must be zero. If, there- 
dei When the Configuration of the Conical Shock Wave Is fore, 
— Known.—Let f, which must be a periodic function of = (6) (13) 
6, take the form wa . 
. : ‘ is the equation of the conical nose, then it is required 
f=e-+ By >oby sin n@é (S‘) I ‘ : . " 
; = that on this surface [Eq. (13) ]: 
where % is independent of 6 and n is odd because the oF oF 
flow phenomenon is symmetrical with respect to the -—— ( )y =f (14) 
: plane 6 = 2/2. Then the description of the flow on Og sin® ¢ \O@ 
> Pane ‘st i¢ fj » r Rae Q avi » rep ¥ , ss 7 
= surface & is defined by Eqs. (9), having — Now in the semiplanes 9 = +7/2, the partial deriva- 
o Eqs. (11) and (12), whenever Yo pre Pi ns “we tive OF /08 is zero and so is ¥, because of the symmetry 
selected (that -" the condition of the a ne "of the motion and of the nose with respect to these 
stream having. been fixed), and, consequent y, & can semiplanes. Hence, for these planes, 
i be taken as the surface S, discussed in (A-1-2), for 
’ ° “ae . ~ 
which the initial data are known. OF /d¢g = 0 (15) 
H | The drawing of the hodographs of the motion in the | nod. 
various meridian planes must be continued up to the 
ted intersection of these planes with the contour surface Upon determination of the value of ¢, by means of 
Lee ‘ . " ° : : ‘ ° 
a of the body; such a trace is determined by the condi- the hodograph in the semiplane 6 = 7/2, such that 
OWS > x ° a . ° . 
Eq. (15) is satisfied, the intersection (generatrix) of the 


a ee 


ues 


tion that upon reaching it the velocity of the fluid 
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nose with this semiplane is defined by the value of ¢ 
so obtained. Meanwhile, the other generatrices, 


which are intersections with the nose of the other semi- 
planes, are obtained more simply by proceeding step 
by step from relation (14). If g = y(6;) is the value 
of ¢ which corresponds to the semiplane 6 = 4,, then 


in the semiplane 6 = 6; + A@ the following is 


(S") 
Og = Fi 


— 0=0; 
¢i4) = Gi + SIN’ G Ad 


(30) 
06 Jp = 05" 


Variation of ¢ and 6, and thus change in the shape of 
It must be 


achieved: 


the shock wave generate shapes of nose. 
noted that through this procedure 0F/0¢ = O must 
again occur, in the semiplane 6 3r/2, for the value 
of ¢ obtained by repeated application of Eq. (14’). 
This requirement will certainly be fulfilled if the de- 
viation of the streamlines occasioned by the shock is 
small enough that the change of entropy from one 
streamline to another may be neglected; this change 
is of the order of magnitude of the product of the square 
of the average deviation by 8;. Under this condition 
the motion can indeed be regarded as irrotational, and 
no contradiction can exist between the assumptions 
and the results to which they have led. 

(A-2-3) Determination of the Conical Field for a Given 
Nose. If the configuration of the tip of the nose is 
given, and considering it to be a right circular cone, 
then the character of the stream at its surface may be as- 
sumed a priori, and this surface of the body may be 
taken as the one possessing the initial data, provided 
that the angle 8; that Vy makes with x is sufficiently 
small in order to make it legitimate to neglect v9°/c* 
in comparison with unity in Eq. (3). Since 0°//06° 
is already of the order of magnitude of 6), the above- 
mentioned assumption leads to the suppression of terms 
of the order of magnitude of 8,°. Indeed it is evident 
that with the above hypotheses it is possible to set 

O°F 


= (=*) + (1/2) sin 2 (2) ‘ 
ii ¢ —— 2) sin 2y | — oe 
; O¢”, o¢ og" 


2 sin? g (F) = 


4 


0 (3’") 
in the immediate vicinity of the body, since v, = 0 on 
it. But expression (3’’) is the differential equation 
that defines the spherical harmonic functions of first 
order. The general integral of this equation which 
satisfies the boundary conditions on the conical nose 


1S: 
F = Ay(e) + Ai(¢) sin 0 16) 


where Ay and A, are a linear combination of the spheri 
cal harmonics of the first and second kind and of their 
associated functions, which need not be specified at 


present. Indeed, with proper selection of the constant 
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coefficients of such a linear combination, it is possible 
to elicit the relations 


dA dg = dA, dg = 0 


which hold for the value of ¢ corresponding to the con 
cal nose, in order that the condition of tangency at the 
body, v, = 0, be fulfilled, while it is also true that 


v7, = An +A, sin 


ve = (A;/sin ¢) cos 6 17 


Now the coefficients Ay, A;, and sin ¢ assume con- 
stant values on the surface of the cones, and, therefore. 
once the values of Ao and A; are fixed, the values 
of the velocities on the conical contour of the body are 
consequently fixed also, and this surface can then hy 
The 
tracing of the hodographs in the various meridian 
planes must be carried out this time until the shock 
The determination of the traces 


taken as the one possessing the initial data. 


wave front is reached. 
of the shock wave in the various meridian planes and 
of the flow conditions upstream of the shock can be ac 
complished by a process of successive approximations, 
which will now be set forth. 

Consider first of all the semiplane @ = 0. 
first approximation, one may let 


In it, asa 


Un = —%MSiNg 
OF | OF f OF 
[4 FS" = - ‘ — = es 
Oe singOodé sing OOo 


(the approximation results because (1/sin ¢)(OF/0# 
is of the order of magnitude of 8), and likewise so is f 


and 


9 


1 — v,? = 1 — %? + %* sin? g = 1 — u%’ cos’ ¢ 


Thus, through means of the equations of the shock 
[Eqs. (9) and (10) ], there arise the expressions: 


— a er 
— sing = SS — Y* COS” 
O¢ k + ] (18) 
F = Up COS ¢ 
Substituting 
% = F/cos¢ (19 


into the first part of relations (18), there is obtained: 


* (tan ¢) al. (1 F*) 20 
(tan ¢ de bal 
And since F and OF/0¢ are known, in consequence of 
the construction of the hodograph, for. any value o! 
¢, Eq. (20) affords the means of determining the value 
of y, satisfying this equation and thus defining the 
generatrix (intersection) of the shock cone in the semi 
plane 6 = 0. The desired value of v for this case 1s 
given directly by Eq. (19). 

Now consider the semiplane 6 = 2/2 (in which f = 
0). Since, in this plane, 


ssible 


coni- 
it the 
‘ 


l7 


con- 
fore, 
alues 
V are 
‘nl be 
The 
idian 
hock 
races 
and 
e€ ac- 


ions, 


10ck 


(1S) 


(19) 
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vz2 = OF /d¢ 


dn = vo (— sin ¢g + sin B, cos ¢); 


and 


1 — v,7 = 1 — %" cos? ¢ — %" sin 2¢ sin B, 


it follows that 


: OF 
sin ¢ + sin By Cos ¢y) ~ ay _ 


fol 


¥ 
=t,. ed ee 
(1 — v%* cos? g — %" sin 2¢ sin B; 
k+1 
and 
F = w (cos ¢ + sin £,; sin ¢) (21) 
From the second of the relations in Eqs. (21) one ob- 
tains 
(22) 


— cos ¢] (1/sin ¢) (22 


sin 8; = [(F/v) 


which, when substituted in the first of the relations (21), 


vields: 
F OF 
w(— sin g + cot g — — cot g cos g) — = 
vo Y 


. 


k— 1 f Bi tad 
b é | l = Vo" cos" a= 29? cos g~ + 2%" CC »s° ¢ = 
ak 9 


Vo 


nae (1 + w* cos? gg — 2cosgm%F) (23) 


Now, since v is known and since the hodograph gives 
the corresponding values of F and OF/0¢ for every ¢, 
Eq. (23) affords the means of obtaining the value of 
y, Satisfying this equation and thus defining the genera- 
trix trace of the shock wave in the plane under consid- 
eration. In the meantime, Eq. (22) permits the calcu- 
lation of the value of {). 

Che generatrices of the wave, or intersections with 
the other semiplanes, are obtained through repeated 
application of Eq. (10), which may be rewritten in the 


form: 
ist hh « 
M%(— sin ¢g sin 6; sin 6 cos ¢) —- = => p 4 
Oy 2 + 1 
(1 — w%* cos? ¢ — v* sin 6, sin @ sin 2¢) (10’) 


Since v and 6; are known and since the hodograph 
furnishes the value of OF /O¢ corresponding to every 
¢g, in the arbitrary general semiplane @ = const., Eq. 
10’) affords the means of obtaining the value of ¢, 
satisfying this equation and thus defining the genera- 
trix for the general case. To check the degree of ap- 
proximation (introduced through the use of the simpli- 
fying hypotheses), it is necessary to verify that the 
first of Eqs. (9) is also satisfied by the values of ¢ which 
are obtained from Eq. (107). 

The above-described procedure presents a first ap- 
proximation for the case where the shock wave and 
the conditions of the flow upstream of the shock are to 


be determined. In general, the accuracy will be ade- 


quate if one stops at this first approximation, provided 
that the magnitude of the angle of attack, 6, is suf- 
ficiently small. If a second approximation is de- 
sired, however, the values obtained thusly through use 
of the first approximation may be substituted into the 
equations of the shock, Eqs. (9) and (10), and the sec 
ond approximation may then be carried out by means 
of a procedure that is quite analogous to the one pur- 
sued in obtaining the somewhat less exact first 
results. 


(B) DETERMINATION OF THE FIELD AROUND THAT 
PORTION OF THE Bopy WHICH LIES BETWEEN THE 
CONICAL NOSE AND THE WING 


(B-0) General Equations 


(B-0-0) Equation of Motion and Associated Charac- 
teristic Surfaces——For the determination of the field 
around that part of the body which follows the nose 
and precedes the wing, the ‘“‘method of characteristics’ 
is applicable. The equation that defines the potential 
of the field corresponding to this region is, using the 
cylindrical coordinates (x, vy, @), 


ra) *Py Pa) *P» 4 10 “Po 4 . re) “Py 
ay . an >= 33 — = 2d) + 
Ox? oy? v? 06? Ox Oy 
] O°?» 1 O°, Vy 
2413 val. a(193 == () 24 
y 00 Ox yoeoy y 
where 
v;" UV, vy” 
ae eee + ne iz =1—-— 
ce co ¢ 
(25) 
U,V¢e v,Ve ve~ 
eS Ss anna @a = 1 -— 
c* c* c* 


First of all, the characteristic surfaces of Eq. (24) are 
sought. To this end, the equation of such surfaces is 
expressed in the form 


W(x, y, 0) = 6 (26) 


where } is a constant. Then IV must consequently 
satisfy the equation :* 


> _ or 6 
- Ox; Ox; 


in which we have set dx; = dx, dx, = dy, and dx; = 
y de. 
The discriminant of Eq. (27) is equal to 


L — [2 + 0,2 + 06°)/c4] 


(for's; & = 1, 3. 3) (27) 


and it is therefore negative everywhere; hence, the 
surfaces defined by Eq. (27) are real. Through solu- 
tion of Eq. (26) for y and thus by transcription of the 
equation of the characteristics into the form 


y = f(x, 6, d) (26’) 


it follows that Eq. (27) may be expressed as: 








$18 JOURNAL OF THE AERONAUTICAL SCIENCES~—JULY, 1949 
cs * 2013. » 223 . 433 + - ——- re! ~ 
Qnjz* — 2aufe + Qn + —J-J¢ — fo + —fe2 =0 
y y r ss 
(27°) 
r 
where /, = Of/Ox and fy = Oof/06. [ 
_ do, 
Now the direction cosines, relative to the straight a we 
dx 


lines x1, X2, X3, of the normal to a general arbitrary one 
of the surfaces [Eq. (27’)] (corresponding to an arbi- 
trary value of the constant 5) and erected at a given 
general point, are: 


—1/V1 + fe? + (fe?/y?) 


fo/yV 1 + fe? + (fo?/y?) 
where the lines x;, x2, and x3; are directed, respectively, 
in this manner: x, along the x-axis, x, along the normal, 
y, to x in the given meridian plane, and x; along the 
normal to the meridian plane. Likewise, the direction 
cosines of the line of action of the fluid velocity in the 
same point and relative to the same lines x, Xe, 3 


are: 
P age 
0.) ¥ Ue" 


9 


+ + 


9 


9 i. 
Vy Ur 
Ve V0, 


ts! 
= 
i) 


Consequently, it is seen from Eq. (27’) that the normal, 
n, to an arbitrary general characteristic surface is in- 
clined to the direction of the velocity, V, by an angle 
the magnitude of which is given as 


— => 


cos' (ns, VY) = = C/V = & om6 
if 6 denotes the Mach angle. That is to say, the charac- 
teristic surfaces associated with Eq. (24) are the en- 
velope surfaces of the Mach cones which propagate, 
from each point, the elemental disturbance created 
there. There are two sets of such characteristic sur- 
faces; one of them (2X) is characterized by the condi- 
tion that the normals to the surfaces make the Mach 


angle 6 with respect to V, while the other one (2’) 


characterized by the condition that the normals to the 


is 


surfaces make the angle —6 with respect to V. 


(B-0-1) Variation of the Velocities on the Characteristi« 
Surfaces.—Consider any characteristic surface what- 
ever of a set (for instance the set =), and designate by 
do, and doz the linear elements of the intersections of 
this surface (denoted by the symbol ~) with an arbi- 
trary general meridian plane, 6 = const., and with a 
plane normal to the x-axis (x = const.), respectively 
(see Fig. 5). 














JL oy 


Fic. 5. Representative trace of a characteristic surface in a 
meridian plane and in a plane perpendicular to the x-axis. 


The velocity changes—dv,, dv,, and d(yv9) 
which occur in traversing the distance do; (the pro- 
jection of which upon x is denoted by dx), are given 


by: 
dv, = (bu + nae (28) 
doy = (Pir + Praf)dx " 
d(yve) = (pis + posfr)dx 


where Py, = O°:/O0x", pro = O°/Ox Oy, and py; = 
0°*b,/O0x08, etc. From Eq. (28) it is evident that 

Pi. = (dv,/dx) — poof, 

bu = (dv,/dx) — f,(dv,/dx) + poof,” 

Psu = (d(yve) dx) — Poafr 
so that now, upon substitution of these expressions into 
Eq. (24), there evolves: 


dv, dv, * 2a13 d(yve) 
ay + (Qufr + 2a12) + 
dx dx y dx 
— . 23 43 ; 
Pr» (Qufr- a aoa. — 2aj2fz) + Pox Z = 2 Te i 
y y 
P33 z 
Je -=0 (29) 
y- y 


Analogously, the velocity changes that take place along 
do» are: 
dv, = (Profo + pis)dé 


dy, = (poof g + po3)d9 


d(yve) = (p33 + Posf9)dé 


where 6 is the angle between the meridian planes within 
the confines of which the said element is contained. 
Whence it follows that 


pos = (dv,/d0) — pofe 
P33 = (d(yve) dé) = fo(dv, d6) +4 Poof o” (28°) 
pis = (dv,/d0) — Priofe 


th 


co 


wh 
Mi é 
bor 
sin 


na 


TO- 
ven 


28) 


nto 


YQ) 


ong 


thin 
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Substitution of these values of p23 and /33 just found into Eq. (29) gives: 


dv, 


s : dv, 2413 d( Ve) 133 d( VUp) 33 « dv, dv, (123 (43: 
mot + (anf, + 2a) — + == ae ft + —¥ (2 _ 97) 4 
wt dx _ dx y dx y? dé y * ’ 0 dé y 


y y" 


* 9 9 
vy hie , . 133 +, =123 » 113 . + 
+ px» (anf + de — 2anf, + —fe’* — lev : ii-) = 0 


bf y y \ 


which, by aid of Eq. (27’), reduces to: 


dv; e dv, 13 d (yup) 33 d (yup) 33 . Qo3 (3 dv, vy, 
+ (—~tal, + Se) —" + 2 Be) fp hoe ee ee Se gee weg s +‘ = 0 (30) 
, dz y = y" y s" jae t. F 


ai 
dx dx y 


This is the relationship that correlates the changes in velocity which take place when proceeding along the two 
line-elements that are the intersections of a given characteristic surface with the two coordinate planes (@ = const., 
x = const.) passing through the same given point. Thus, Eq. (30) is the relation that corresponds to the one in 
two-dimensional motion that states the change in velocity along a characteristic element and is equivalent to that 
case in all respects with the mere substitution of the system of two orthogonal elements, shown in Fig. 5, for the 
single linear element with which one is concerned in the two-dimensional case. The element contained in the 
plane x = const. corresponds to the term indicated in Eq. (30) by means of the symbol [],,, which represents 
therein the influence upon the fluid motion of the fact that the flow is three-dimensional. When this relation is 
applied to the determination of the potential field, it is readily seen that this element does not introduce any 
greater complication than in the case of two-dimensional flow, at least from the standpoint of the concepts in- 
volved, because this element is upstream (with respect to the direction of motion) of the associated element belong- 
ing to the other orthogonal plane and therefore the conditions of flow on it are known independently of the condi- 
tions existing on the associated element. 
If, now, Eqs. (28’) are substituted into Eq. (24), then by taking into account Eqs. (28), it follows that 


dv, dv, e 2ai3 * 2a3 dv, 2d23 (133 dv, (133 d( VU) Vy 
1 + 2ayo — Ayn fr — he )| + | —— — ——Jj + 8) See = Q (30’) 
[ " dx dx ( rr sa y x y do y y? °} do rr @s 7 


So that from comparison of Eqs. (30’) and (30) it results that the following must be true: 


| diy ; to) n E in diy .| 7 (31 
dx fo i i. dé d@ Je ee ‘a i 


his constitutes another relation between the velocity changes under consideration which take place along the two 
associated elements belonging to the two orthogonal planes, and it is substantially a necessary consequence of 
the irrotationality of the motion. 

(B-0-2) Fundamental Equations of the Characteristic Lines.—Relations analogous to those of Eqs. (30) and (31) 
hold for the changes of velocity which occur along the elements do,’ and do’, which are the intersections of a char- 
acteristic surface of the other set (’) with the planes @ = const. and x = const. In view of this, we now let f, = 
uw at an arbitrary general point and denote by \; and z the two values of /, which are obtained from Eq. (27’) 
u is used therein (they are necessarily real if u is sufficiently small, and this certainly happens in the 
Then the following fundamental equations hold for the characteristic lines (which is 


when fg = 
case under consideration). 
the name given to the intersections of the characteristic surfaces of the two sets with the planes 6 = const. and x = 


const.) : 
(A) For the characteristic lines 0; and o2 belonging to the set of surfaces 2: 


4 = ¢; fA = Ay 


| (= + | 4 | 2% dv, 4 (= “) dv, 4 33 fone] 4 Vy 0 
2 2 — 3: = 
_ dx * GF do. y do y aed y*/} dé y> dd 1, Y (32) 
dv, Hon) | |e | 
—- _- — : N = 0 
| a” ah’ ae "Si 


where the third of these expressions is obtained from Eq. (30’) through utilization of the relation between the roots 
\i and d, of Eq. (27’). Now denote by 72 the tangent to o;’ at the point P, which is considered to be in the neigh- 


borhood of the line-elements do; and do,’ of the two characteristic lines lying in the plane @ = const.; then dy = 
Sil (72, x)/cos (72, x). And upon letting 


dv; = dv, cos (T2, x) + dv, sin (r2, x) 
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further reflection shows that the third of Eqs. (32) yields: 
Qy3 dv Ao; dv, 33 d(yve) d?; v 
[2 oe (2 3 Qs “) dVy ea = i) cos (rs, x) + (au —*) + — [cos (rs, x)] = 0 (32) 
y dé y y?/ dé Ff @ i. dx/., Y 


It is easy to see that (dv,;/dx),, defines the variation, along the element do,, of the component of velocity that is 


yarallel to the other element do,’ passing through the same point P. 
} Pp g 8 I 


Written in the form of Eq. (32’), this rela- 


tion gives an even more precise appearance to the analogy between it and the relationship that holds for the case of 


two-dimensional flow, in which one has simply: (dy;/dx),, = 0. 
(B) For the characteristic lines 0,’ and 2’ belonging to the set of surfaces 3’: 


to = B; 


dv; “| | =" dv, 2d23 m ) dv, . a33 d( me Uy 
be ee + (= - oe +S | + Fao 
[ ( ili dx By y dé y " y’/ dé of y> db Ja y (33) 


| dv, ty | 4 
dx” dz i» 


E 4) or ait 
dé * 0 Le 


Again, denoting by 7, the tangent to o; at the point P considered above, and then letting 


A, = sin (7, x)/cos (7, x) 


and 


dv, = dv, cos (1, x) + dv, sin (7, x) 


it follows that the third of Eqs. (33) will yield: 


2a33 dv A»: dv, 33 d( Vv») dv v 
| =" —o + (2% — = “) y .Y = BAd) cos (71, x) +. (a: “n + v cos (71, x) = () 
y dé y y’/ dé y? dO Jy dx/, YY 


(33’) 


v1 


It is apparent that (dv,/dx),, defines the variation, along the element do,’ of the component of velocity that is par- 


allel to the other element do, through P. 


Assume that the two line-elements P,P and P;’P of 
the characteristic lines o; and o,’, respectively, that 
pass through the point P, are given and represent on 
the hodograph plane (v;, v,), corresponding to the merid- 
ian plane under consideration, the two vectors defin- 
ing the velocity components at P; and at P,’ which lie 
in this same plane. Then Eqs. (32’) and (33’) make it 
possible to determine immediately the analogous veloc- 
ity component at P, provided the flow conditions are 
known on the characteristic surfaces © and 2’ to which 
the points P; and P,’ belong (or at least known in the 
neighborhoods of the surfaces that contain P; and P,’). 
Under these conditions it is possible to calculate the 
values of dy; and dv, from Eqs. (32’) and (33’), with the 
aid of which the graphical construction indicated in 
Fig. 6 is derived. 

This construction is an obvious extension, to the 
problem of three-dimensional motion under considera- 
tion, of the procedure used for two-dimensional motion 
and suggested by R. Sauer for the case of axially sym- 
metric flow. 


(B-1) Field Corresponding to the Expansion Resulting from 
the Change in Direction of the Meridian Line 


(B-1-2) Delineation of the Characteristic Surface at 
Which the Expansion Begins.—The application of Eqs. 
(32) and (33) to the determination of the field of flow 


around that portion of the body that lies between the 
conical nose and the wing, when the field around the 
nose is known, is carried out by means of a step-by- 
step procedure that is entirely analogous to the one 
that is followed for the case of two-dimensional flow 
or for the case of axially symmetric motion. 

Let AA be the circumference of the body (for a sec- 
tion normal to the body axis) where it adjoins the coni- 
cal nose. Then the perturbation due to the change in 
direction of the meridian line begins to be propagated 
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(a) (b) 
Fics. 6a-6b. Graphical construction for determining the 
velocity at the intersection of two characteristic surfaces when 
the velocity is known on them nearby. 
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within the fluid starting from a characteristic surface 
passing through AA. This surface is denoted by 
>, and it is determined on the basis of the results of 
the investigation discussed in Sections (A) and (B-0). 
Let the equation of this surface have the form 


Y = Va + ao (0)(x — x4) 


in the neighborhood of this parallel AA (the above- 
mentioned circumference); then (u),-2, = (fe)zar4 = 0 
for any 8, and \y = 5 = ap (8). 

Upon substitution of these values into Eq. (27’), it is 
possible to calculate at once the values of a)“ = dA, 
which correspond to the various meridian planes, in 
consequence of their being solutions of the second de- 
gree equation @y,A? — 2ay~4 + dx = 0. The coef- 
ficients of A in this equation are known because the 
velocities are known. Let 6 denote the angle of inclin- 
ation, with respect to x, of the component of velocity 
lying in the meridian plane, and let 6; be the local Mach 
angle corresponding to this component (or sin 6; = 
¢ V 0,2 + v,*); then it may be stated exactly that 
\, = tan (B + 6,) and A, = tan (8 — 6,). Thus, it is 
possible to draw the tangents to the lines o; (which lie 
in the various meridian planes) for the point at which 
x = x, and y = y,. These tangents determine the 
tangent cone at AA for the surface © under considera- 
tion. Upon cutting this cone with the plane x = 
x, + Ax, therefore, the line with equation 


Xa t+ Ax = Xa, 
Ya + ao Ax = yu, 


x 
Vv 


is obtained for a sufficiently small value of Ax. In 
the neighborhood of this line, let the equation of the 
same surface be given in the form 


yY = Va, + a(O)(X — Xa,) 
and therefore it follows that 


(MH)z=za, = (fo)eaxa, = Go Ax 


= a,” (8) 


A _ (r)x =A, 


Consequently, when the values of /» and f, just deter- 
mined are used in conjunction with Eq. (27’), it allows 
the determination of a;(6) for the surface ©, and in 
this manner the whole shock front of = is determined 
step by step. 

The actual meridian line of the body is then replaced 
by a line of jointed segments inscribed in the original 
contour in such a way that the change in direction at 
each corner of this line is small enough to be propagated 
by only one characteristic surface of the set 2. Letting 
~" be the surface of this set which propagates the per- 
turbation originating at AA; then the equation of 
="), in the neighborhood of this parallel, is in the 
form: 


Y= Va aa a,’ (A)(x — X4) 
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whence again, (u),-2, = (fe)rmza = 0, for 6, and \; = 


a = a, (8). 

In an arbitrary general meridian plane, consider the 
change of velocity along an element Ag,’ of characteris- 
tic line (belonging to set =’) that passes through a 
general point Po of o;, which is the intersection of 
the meridian plane in question with the surface =. 

This variation of velocity satisfies the third and fourth 
of Eqs. (33). If now Av,, Av,, and A(yve) denote the 
changes in the respective quantities between P) and 
the point P,;"? where the element Ao,’ intersects o;‘” 
(Fig. 7), it follows that (here, for sake of simplicity of 
notation, the terms in the brackets [],,,, given in Eq. 
(33), have not been repeated): 


ay( Av, 4. Ai Avy) oy" + Ax [oy + Ax(v,/y) vid of (34) 
A(yv6)e," ed Ax[ Joy’ = 0 


in which ay, v,/y, and the coefficients of Ax must be 
calculated at Po and for the element of the character- 
istic line 2’ which is normal to the meridian plane con- 
taining Po and P,™. 

Now if we pass to the limit, as Ax — 0, and therefore 
as Po and P, are made to tend to A along their 
respective characteristic lines (Fig. 7), Eq. (34) reduces 

“~ 
( Ave )iim a _»,\ 


to 
li (=) 
im = 4 oar 
Ax=0 Av, oi’ ( Av, )tim ( (34 ) 
A(yve)o, = 0 





The second relation of Eqs. (34’) states that component 
ve does not undergo any discontinuity in passing from 
> to 2 at the corner AA. The first of Eqs. (34’) 
expresses the important property that, in every me- 
ridian plane, the changes in the velocity components 
Av, and Av, at the parallel AA follow the same laws 
that hold for two-dimensional motion,! since nh = 0, 
and therefore \; is simply determined by 


QyA?2 — 2ZaprA + doe = 0 


For the purpose of calculating these changes in veloc- 
ity propagated by the wave >‘” and corresponding to 
a given change in direction Af, one may use the formu- 
las for two-dimensional motion :° 


ss 
Cx 


( Av) sim — a3 ( -} a) 


f' » 
Vv M,? - 1 
vy ) 
= 
V My = l 
where v,“ and v, are the velocity components at A 


on >“ and Mp is the local Mach Number at the same 
points, relative to the velocity component in the merid- 


( Avy)uim = AB ( v, 


ian plane only 


{an = cet een 


The first of Eqs. (34’), on the other hand, permits the 
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Fics. 7a—-7d. Extension of the hodograph from meridian plane 
to meridian plane for the expansion occurring at the base of the 
nose cone. 
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(e) 


Fic. 7e. Determination of the velocity at ?," and Q,\ by 
use of computed components of the increments that are directed 
along the characteristic lines o;’ and o;""’. 


calculation of \;, and therefore the determination of 
the tangents to the characteristic lines o;‘” in the vari- 
ous meridian planes. 

(B-1-2) Determination of the Characteristic Surface 
Marking the End of the First Expansion, As Well As 
Description of the Velocities Obtained As a Result of 
the Expansion.—Now consider any point P,‘" on the 
characteristic line 0," just determined, for the plane 
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6 = 7/2; it also lies on a characteristic line o,’, which 
passes through P,‘” in the same plane, and it is close 
enough to A (Fig. 7) so that the are Po P,"” of o;' is 
indistinguishable from its tangent (which is known) at 
P,° This Section 
(B-0-2), and the velocity at P,“ 
through application of the graphical construction that 


is the situation considered in 


can be obtained 
was described there (Fig. 7b). In consequence of the 
above, it is possible to find the tangent to o,°” at P, 
by aid of Eq. (27’), with fy = 0 therein, and it may now 
be drawn in. 

Then consider, in the meridian plane 06 = (7/2) 4 
Ad, the intersection, Q,"'’, of the tangent to o;") at A 
(already drawn in, following what was explained above 
with the plane x = const. = xp,‘”, to which the point 
P, belongs (Fig. 7c), and let Qo“ be the point of the 
line o;°", which belongs to the same characteristic line 
a,’ that passes through Q,"”. 
make it possible to determine the velocity 


sy? 


In turn, the Eqs. (32 
and (337) 
component at Q,"', lying in the plane @ 
Ad, while the fourth of Eqs. (32) [or the fourth of Eqs. 
33)] allows the calculation of the difference between 
the values of (yve) at the points Q;‘" and A [or at the 


points Q,"? and Qo] and therefore the determination 
1 


(r/? + 


of the value of the same component v¢ at Q,‘” (see Fig. 
7d; in the same Figs., 7a, 7b, and 7c, it is shown how the 
values of (dv,/d@),,, etc., corresponding to the linear 
elements of the characteristic oo’, are obtained). 
Writing the equation of Y‘' in the form 


y= va + [xe xa lay (6) + [x — xpi’ Jae! (0 
in the neighborhood of the line that is the locus of 
points P;“”, Q,”, ete., it follows that 

u)> = (foenrp(? = (xp — xa)" 
which is known as a result of the preceding discussion 
Then, knowing the velocity components at Q,"”, it ts 
possible to calculate the value of Ay = a2"?[(/2) + 
Aé| at QQ)” by means of Eq. (27’), and thus it is pos 
sible to draw in the tangent to the characteristic line 
a," at the same point Q,'”. 

Having calculated the velocity at the points P, 
Q,'", ete., and having determined the tangents to 
o;‘") at these points in all the meridian planes, one 1s 
again in a position to calculate the velocity at another 
P,". 
procedure just explained makes it possible to obtain the 
whole S“ and [through depiction in the hodographs 
related to the various 0; "| the velocity at all its points 


In this manner, the repeated application of the 


The following is a practical application of the method 
outlined above. 

Table 1 lists the values of v,, v,, ve, dv,/d6, dv,/ dé, 
and dve/d@ at such points as Py, Qo", etc., belonging 
to the various meridian semiplanes. Table 2 presents 
the same information for the points A belonging to 
o;"'’ for the various meridian semiplanes. 
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TABLE 1 


vhich 
lines Values of the Velocity Components and Their Derivatives at Such Points as Po), Qo), etc., Belonging to the Various Meridian 
. — Semiplanes 
itis A x ; dv; diy dup 
n) at ' ad "y ” d0 dé do 
ection 90 0.6385 (0.09593 0 0 0 —0.04799 
ained 70 0.6380 0.09480 0.01679 0.00252 0.00685 — 0.04503 
50 0.6367 0.09154 0.03154 0.00470 0.01243 — (0.03660 
that 30° 0.6346 0.08658 0.04254 0.00634 0.01722 —().02407 
f the 20° 0.6334 0.08360 0.04621 0.00696 0.01922 — 0.01657 
Po 10° 0.6321 0.08034 0.04852 0.00709 0.02071 — 0.00844 
0° 0.6309 0.07692 0.049388 0.00701 0.02116 — 0.00016 


7 now 


I 
i 
TABLE 2 
at A Values of the Velocity Components and Their Derivatives at the Points A Belonging to o;“) for the Various Meridian Semiplane 
' 


bove ‘ a dv; dv, dvg 
; u Y% "0 a 
point : ; dé dé dé 
f the 90 0.6710 0.07447 0 0 0 — 0.06266 
70° 0.6703 0.07440 ().02143 0.00351 0.00039 —0. 05888 
¢ line 50° 0.6685 0.07420 0.04028 0.00659 0.00073 —0.04800 
(39’ 30 0.6657 0.07389 0.05426 0. 00896 0.00102 — 0.03133 
- 20° 0.6641 0.07370 0). 05888 0.00980 0.00112 —0.02143 
locity 10 0.6623 0.07350 0.06171 ().01017 0.00113 —0.01088 
9) 4 0° 0.6605 0.07331 0.06266 0.01011 0.00112 0 
| Eqs. 
[weet ° 
it the On the semiplane 6 = 90°, the following values may istic o;"° through P™. From the end point of this vector 
ation be calculated at Po: draw a line perpendicular to it [and therefore perpendic- 
e Fig. 2 oe ular to o;°]. These two perpendicular lines intersect 
‘a an = —2.496; a3; = 1.000; Ax/y = 0.06699 1] per] 
roclaons. Pit, _ ata point O,. The vector O,O, is the change of veloc- 
linear assuming for y the average — between the ordi- ity in passing from A (of 0") to P,. From Fig. 7e 
bees . : . (0). (1). ane oN se EIT = on : 
nates of the points Po and Py 5 CoS (71, %) = 0.7607; it is seen that this velocity change has the components 
thus the equation that determines Av, is Av, = —0.00282 and Av, = —0.01181; thus the veloc- 
—().04799) (0.06699) (0.7607) + (—2.496) Av, + ity components at P; are v, = 0.6710 — 0.0028 = 
9 . ; } 


(0.09593) (0.06699) (0.7607) = 0 90.6682 and v, = 0.07447 — 0.01181 = 0.06266. The 
same procedure can be employed for 6 = 0°. At the 


‘us of thence » = OR . . ° ° 
whence Av, +0.00098. point Q)“ in this plane, the following values may be 


At the point A of the same semiplane belonging to etalon 
o\\", the following values may be obtained: 
en Qy = — 2.352: a3 = — (0.2624; ao = — 0.03199 
) dy = —3.136; a33 = 1.000; Ax/y = 0.5756 2 
ssion. ff , ; ad33 = 0.9795; = —0.01338; cos (7, x) = 0.7678 
tj along the segment AP,” taking the average of the y 
1 S 
/ ordinates at A and P,‘"]; cos (72, x) = 0.9210; thus Ax x 
2) + de | “ — = 0.06214; the term [],,, = —0.00551 
the equation determining Av; becomes y 
S pos - 
ic line 3.136 Ave + (—0.06266) (0.5756) (0.9210) a so that Av, = 0.00145; at A: 
. (0.07447) (0.5756) (0.9210) = 0 7 aa ae 
P : a = —2.907; ayy = — 0.3706; ade = — (0.041135 
rte 
whence, Av; = .00200. ; uM 1 dy . 
ae: ence, Ave = +0.0020( a3 = 0.9648; ~ = 0; 7 = —0.00105 
. aN] dt 
one 1S Indicating the components of velocity at A with the ae ° 
1other subscript “1’’ and those at Py“ by the subscript “2”, - = 0.5490; cos (72, x) = 0.9146 
, > . Sal ; . . . 5 
a the consider the vector a whose components in the meridian 
in the plane are v7}, — 2%, = 0.6710 — 0.6385 = 0.03249; whence Ax, = 0.01128. 
rT: Ss E —- —_ - . e ° . ° ° 
raph : Vy — Vey 0.07447 — 0.09593 = —0.02146, respec- Repeating the graphical construction indicated 
points: F tively. This vector is represented in Fig. 7e by the above as shown in Fig. 7e, the following values are ob- 
> h d ' — _ “ ens ' ° 
cum directed segment OO,. Add to the vector @ a vector of tained: 
j the same magnitude as Av; and parallel to the charac- Av, = —0.00051; Av, = —0.02906 
lv ,/ a0, teristic line o,’ (through A) and draw through the end ' : 
ying ° ° = ¢ 1) , —_ C 
miging point of the vector sum the normal to the vector drawn thus at Qi“ (for @ = 0°), 
se parallel to o,’.. Analogously, draw through the origin v, = 0.6605 — 0.00051 = 0.6600 
ing to 


0 of @a vector equal to Av 


né 


ind parallel to the character- vy, = 0.07331 — 0.02906 = 0.04425 
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Fics. 8a-8d. Extension of the hodograph from meridian plane 
to meridian plane for the development along the straight-line 
portion of the first frustum. 


In order to calculate vg at the point under considera- 
tion the fourth of Eqs. (33) may be used. This can be 


written in the form 


Ax | dv, dv, ) be 
A» =i% = 20,” 
y E + 5 (3 , | v y 


dy/dx along o;’ being equal to A, and Avg being the dif- 
ference between ve at Q,"? and Q)'””. 


Ave = 


From the above values 


Ave = 0.06214 [0.00701 —0.4830 (0.02116 —0.04938) | - 
(0.04425 0.07692) (—0.01338) = 0.00085 


that is, 


ve = 0.04938 + 0.00085 = 0.05023 


(B-2) Field Corresponding to the Development Along the 

Straight-Line Portions of the Meridian Line 

Let BB be the parallel of the body that contains the 
vertices B of the line of jointed segments which is in- 
scribed in the meridian line and which is assumed to be 
close enough to parallel AA that the portion of char- 
acteristic line o;’ contained between B and 5‘ can be 
replaced by its tangent, at the point of intersection of 
o;’ with © in every meridian plane. In the meridian 
plane @ = 7/2, if P,"” is this point of intersection, the 
changes in velocity along P,;B satisfy Eq. (33’) 
Therefore, since the velocities on =“) are known in con- 
sequence of what was found in Sections (B-1-1) and 
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(B-1-2), this equation affords at once the means of cal- 
culating the component Av, [in the direction of o,! 
through P,] of the vector difference between the 
velocities at B and at P,. 

On the other hand, the velocity at B has the imposed 
direction AB, and therefore the construction shown in 
Fig. 8 is evolved for the calculation of the velocity at B 
in the hodograph which is related to the plane 6 = 
1/2. Thus the tangent to o, through B can also be 
traced in at once by use of Eq. (27’), because fs == f 
for 9 = 2/2 and because the velocity components at B 
are now known. 


In the meridian plane 6 = (w/2) + Ad, for the sake 
of simplicity of notation, again let Q,"? be the point 
where the line o;’ (passing through the point B in the 
plane under consideration) intersects the plane >“ 
Then a construction like that shown in Fig. Sb allows, 
in view of relation (33’), the determination of the veloc- 
ity components, in the plane under consideration, at 
the point B. The construction is carried out in the 
hodograph that is related to this same plane, the 
velocity at Q; being*known in consequence of the 
(See Figs. Sc 


derivation presented in (B-1-2). and 


Xd. ) 
Further, because of the irrotationality of the motion, 


it is necessary that 


A(yve) V; = _. ae ap 
= cos Bp + sin 8 oh 


As eli] 


holds, where 8’ is the angle at which AB is inclined to 
x and where A(yvg) is the difference between the values 
of (yve) at B and at A and As is taken as the length of 
AB. Since, now, 

Ov, (zp) 6 = (w/2) +40 (Yzp)o=n/y 


o#4 A@é 


and Ov,/06 are known, the calculation of (yvg) at B, in 
the plane 6 = (4/2) + A6, may be carried out by means 
of Eq. (36). Thus, in turn, the slope of the tangent to 
the o,° at B can be calculated by means of Eq. (27' 
in which (fo) must always be set equal to zero; this 
tangent can then be drawn in in the plane under con 
sideration. Proceeding in a like manner from plane 
to plane and proceeding from point to point in the 
same plane, it is possible to complete the tracing of the 
+“) which -passes through the parallel BB and likewise 
the construction of the hodographs relative to its 
characteristic lines in the planes 6 = const. 

The practical application of the method outlined 
above is made in a manner analogous to that indicated 
in Section (B-1-2); therefore, this application will not 
be repeated either in this section or in the following 


ones. 
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Fic. 9. Representative trace of a characteristic surface in planes 


which are normal to the Y-axis and the x-axis. 


(C) DETERMINATION OF THE FIELD FOR THE REGION 
OF FLOW IN WHICH THE WING Is SITUATED 


C-0) General Equations 


C-0-0) Potential Motion in Cartesian Coordinates 
and the Corresponding Equation of the Characteristic 
Surfaces. Let the potential in this region be #3V;. 
For the determination of the motion in the neighbor- 
hood of the wing it is convenient to make use of the 
Cartesian coordinates (x, Y, Z) and to take the plane of 
the wing as belonging to the coordinate plane (x, Y). 


BETWEEN 


C 
~] 
or 


WING AND BODY 


et 


For a study of the motion in proximity to the body, 
however, it is still convenient to use the cylindrical 
coordinates (x, y, ). 

By use of the Cartesian coordinates (x, Y, Z), the 
differential equation that defines ®; is: 





9) "Ps o *h; Oo “Ps 
A "ea + 2A 2 weal Ax» — - 
" Ox? ” dx 0] + oy? 
0°d; 0°d; Od, 
2A13 —— +24, —— + Ap— =0 (37) 
"ta "trae Se _ 
wherein 
An = 1 -™ v,"/c* = a3 Ap = —v,0y/c* 
Aig = —v,0;/c?; Ax = —vytz/c? 
An = 1 — (vy?/c? ; Ay = 1 — (vz?/c*) 


Analogously to what was demonstrated in (B-0), it can 
be shown that Eq. (37) admits two sets of real charac- 
teristic surfaces. Letting the equations of these sur- 
faces be of the form 

x = f(Y, Z, b) (38) 
where + is again an arbitrary constant, the functions f 
must satisfy the equation: 


Ay — 2Anfy — 2Anfz + Anfy? + 2Anfyfz + 


A3sfz" = ( (39) 


from which it follows that the characteristic surfaces 
may be interpreted to be the envelope of the Mach 
cones that propagate the elemental perturbations, as 
has already been observed in (B-0). 


C-0-1) Variation of the Velocities Along the Characteristic Surfaces.—For any characteristic surface of a set (still 
denoted as S), let do; and do» again be the linear elements of the intersections of this surface with the planes Y = 


const. and Z = const., respectively (Fig. 9). 


As in (B-0-1) and (B-0-2), it is found that the increments of velocity, 


dv,, dvy, and dzvz, taking place along these linear elements, satisfy the equations: 


dv ; dv, lv, lv; ' .. dv, 
E E it. (dead & BAW) oe DAs - | + | As t + i{~bo Oha + Sha: 4 =@ (40) 


dY % d "49 


which is analogous to Eq. (30"), and 


diy . . dv, 
A »” (2A ” Eacegs 2Aos/ = A ofy) 2/ 


which is analogous to Eq. (30). 


dZ dZ 


lvy = * loz . 
fe = + (2Ay — Asfz) - +A | = 0 (40’) 


~ i2 dZ 3 dZ 


Here, dZ is the projection of the element do, upon the x-axis and d¥ is the projection of do, upon the Y-axis. 


From comparison of Eqs. (40) and (40’), it follows that 


| . dv, dtvz 4 EB 4} = | — 
ay a7. tae WA 


which is analogous to Eq. (31). 


(C-0-2) Fundamental Equations for the Characteristic Lines. 


(41) 


Relations that are analogous to those of Eqs. 


(40), (40’), and (41) exist for the changes in velocity along the elements do,’ and do»’, which are the intersections 
of the planes Y = const. and Z = const. with the characteristic surfaces of the other set (’). 
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Let fy = uw, and let \, and dz denote the two roots of Eq. (39) when fy is set equal to yu in it (these roots are neces- 
sarily real if u is small enough, as certainly happens in the problem under consideration); then the following equa- 
tions hold for the characteristic lines o; and o2, which are the intersections of the characteristic surface of the set 


const. and Z = const.: 


Y with the planes Y 


diz 
i, 
| - Ee 


fy = 


J 


a) i‘ | { dvy 
G27 L., el) 


ao 


| . dv, 4 . | 4 
~ aY dY4,, dZ 


For the characteristic lines 0,’ and a2’ belonging to the set of surfaces ¥ 


ty = #; 


' ~ ‘ [ dev 
eri. ees 
| dv, 

— ho _ — 
dJ 


Denote by 72 the tangent to the o;’ at an arbitrary 


+ { 


dvz 
A33 ( 
| ; dZ r 


general point P on it, and let 
dv, = dv, cos (2, x) + daz sin (72, x) 


then since A» = cos (Ts, x)/sin (72, x), the third relation 


of Eqs. (42) produces 


[A33(dv,/dZ) ],, + []o, sin (72, x) = O (42’) 


and, as in Section (B-0-2), it is evident that (dv;/dZ),, 
defines the change along the linear element do; of the 
velocity component parallel to the other element do,’, 
which passes through the same point P. 

Likewise, denote by 7; the tangent to the o, at P, and 
let 


dv, = dv, cos (m1, x) + duz sin (7, x) 


then, since 


A = COS (7, x)/sin (74, x) 
the third of Eqs. (43) gives 


[A33(dv,/dZ) Jo + [Jo sin (m,x) =O  (48’) 
where (dv,/dZ),,, defines the change along the linear 
element do,’ of the velocity component parallel to the 
other element do; which passes through the same point 
r> 

Therefore the construction described in (B-0-2) 
again valid for the determination of the velocity com 
ponent at P and which lies in the plane Y const., 
provided the two segments P,P and P,’P of the charac- 
teristic lines o; and o2’, passing through P, are given and 
provided further that the conditions of flow in the 
neighborhood of the points P; and P,’ on the charac- 
teristic surfaces © and S’, to which they belong, are 


iS 


known. 
It must also be noted that, if there were no interfer- 
ence effect of the body upon the wing, then the term 


ob ( — A sou + 2A 12) 


dvz 
i). 


by; he =X, 


dv, dvz 
2 24,, 
dY ba . 3 | 


dv, (42) 
Rees = 0 
az \,. 
y’, one has, analogously, that 
{ + 2A,) Ws +. 24 dez 0 
é1 90M -/112 dY <4 ay bs 13 


E dv, | 
+ u = ( 
je * aL. 


within the empty bracket in Eqs. (42’) and (43’) would 
be zero, and these Eqs. (42’) and (43’) would reduce 
to the relations dv;/dZ = 0 and dv,/dZ = 0, which ar 
the equations of two-dimensional flow. The construc- 
tion just described would then reduce to the well- 
known construction based on the reciprocal relationship 
between the bounding net of the Mach lines in the 
plane of motion and the limiting (principal) net of the 
analogous characteristic lines in the hodograph plane.’ 
These same terms, which go into the empty brackets, 
therefore describe the effect of interference on the local 
changes of velocity at the various points of the velocity 
field. 


C-1) Description of the Field Existing on the Top Surface 


of the Wing i 

(C-1-1) Expansion Field Around the Leading Edge 
The application of the formulas (42), (43), (42’), and 
(43’) to the determination of the field of flow on the 
upper surface of the wing is performed by means of a 
procedure that is analogous to the one expounded im 
Section (B). The wing is taken to be a flat plate that [ 


lies in the plane Z = 0. Without loss of generality i 
application of the procedure, it is permissible to assume 
that the fluid stream suffers an expansion on the upper 
surface of the wing; indeed, it will become clear as 2 
result of the considerations to follow what the procedure 
should be in case a compression rather than an expat 


sion occurs, or in case the wing, instead of being in the 
form of a flat plate, has any contour whatsoever (pro 
vided that it is compatible with a stream that is every- 
where supersonic). 

Also, for the sake of simplicity in explanation of the | 
procedure and in the related formulas, it is herewith 
assumed that the body is cylindrical, at least over that 
portion of the body which is affected by the wing; it's 
assumed that the radius of this cylinder is R, and the 
plane Y = Rcan be regarded as the plane of symmetr) ; 
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INTERFERENCE BE 
as far as the study of the flow over the wing is concerned. 
This hypothesis makes it possible to determine the flow 
on the wing and body separately, for each of which it is 
convenient to apply a different reference system, and 
yet it is not likely to alter the results of the investiga- 
tion in any appreciable way. 

C-1-1-1) Delimitation of the characteristic surface at 
which the expansion begins: The perturbation produced 
by the wing is propagated into the stream flowing over 
the upper surface beginning with the characteristic sur- 
face arising at the leading edge of the said wing. This 
surface is still indicated here as =, since there is no 
danger of ambiguity; it is possible to obtain its de- 
scription at once, in view of the investigation presented 
in (B). 

If the equation of =“ is written in the form 


x = ZA)(Y) (44) 
for the neighborhood in proximity to the leading edge, 
then it follows that 


(fz)zeo = Ao(¥) =»; (fy)z-0 = 0 


and if these relationships are substituted into Eq. (39) 


there is obtained: 
Ay ns 2A 3X + A3z3\" = O 


Thus it is possible to arrive at the value of \, through 
the utilization in A, of the known values of the veloc- 
ity components on the edge of the wing and upstream 
of the surfaces >. 

It is deduced that A; = cot (8 + 6,), where B now 
denotes the angle of inclination, with respect to the x- 
axis, of the component of velocity in the plane VY = 
const. and 6, stands for the local Mach angle corre 
sponding to this component 


/ 


(sin 6; = c/Vv,? + 7,”) 


In the various planes ]° = const., it is thus possible 
to draw in the tangents to the lines o,;, which are the 
intersections of the surface =" with these planes and 
which constitute a ruled surface that is tangent to the 
characteristic surface at the edge of the wing. 
rhen pass the plane Z = A» through this ruled surface, 


ind the line with equation 


Ww (O) 
—~ 


AnAo { VY) = My 
Ao 


J x= 
\Z = 


is produced. 


Let the equation of = in the neighborhood of this 
line be expressed in the form 


x1 + (Z — Ay) Ay (YY) (45 


x= 


then it follows that 


Uz)z “4, = A, ( ‘é = Ny 
AnpAy ( VY) 


(fy)z 2A = 


TW 
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When the values of fe and fe. just determined are sub- 
stituted into Eq. (39), it affords the means of getting 
A, = A, (Y) and therefore of continuing the tracing of 
the o,;“° in the planes 1 = const. Thus the entire 
surface © is determined by proceeding step by step. 

(C-1-1-2) Determination of the characteristic surface 
constituting the boundary between the first and second of 
the n expansions used to align the flow with the wing sur- 
face at the leading edge: The total change of direction, 
B(Y), that the streamlines must undergo at the leading 
edge, in every plane Y = const., in order to align them- 
selves with the plane of the wing, is assumed to be pro- 
pagated through means of n characteristic surfaces of 
the set =, in such a way that each one of these surfaces 
propagates the change AS(Y) = B(Y)/n. 

The equation of =” in the vicinity of x = Oand Z = 
0 may be written in the form x = ZA,"(Y), with the 
consequence that 


(fz)z=0 =) = A,')(Y) 
(fyz-0o =u = 0 
In the arbitrary general plane }’ = const., the velocity 


changes, along an element Ao,’ of the characteristic 
line (belonging to the set =’) passing through an arbi- 
trary general point Po(°) of o,°° just obtained above, 
satisfy the third and fourth relations of Eqs. (43). 
Consequently, if Av,, Avy, and Av, denote the changes 
in the components of velocity between Py) and Py, 
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Fic. 10. Relation between hodograph and characteristic sur- 
faces for velocity changes through shock at wing leading edge. 
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then we have 
[A 33( Avz , Ar’, lo,’ “V AZ| los! 0 1G 


[ Avy |," AZ[ J. = 0 


Of course Po is the intersection of the straight-line 
element Ao,’ with the line o;“ (Fig. 10). 

Now, therefore, if we pass to the limit, as AZ — 0, 
and thus for Po and Po‘” approaching the leading 
edge of the wing along their respective characteristic 


lines, Eqs. (46) reduce to 


Avz \ Avz)tim 
lim = 
o \Az,/ az 0 Av) um 


AZ 
Avy t 9 


rhe interpretation of Eqs. (47) is entirely identical 
to that pointed out in Section (B-1-1)—namely, that 
the velocity component vy does not suffer any discon- 
tinuity across D at the leading edge of the wing. Mean- 
while, atx = 0 and Z = 0 the changes in the compo- 
nents v, and vz, lying in any plane Y = const., follow 
the same laws that are valid for two-dimensional mo 


( / “: + vz ») 
nN V M,? -— 1 
s } (4S) 
o,,(9) 
B(Y ) (-s, a 4 Uz ) 
l 


n VM? — 


tion. Therefore there is evolved: 


( Av) tim = 
{ Avz ) tim 


where v, 


at the leading edge of the wing and ahead of =“, while 


and vz are now the velocity components 
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Mo is the Mach Number at the same locations and js 


defined as 
M, = VIR FS Pe 


The value of \; may then be computed at once by 
means of Eq. (39) [or by Eq. (47)], and therefore the 
tangents to the characteristic lines 0; belonging to 
=“ and lying in the various planes Y = const. may be 
determined. The evaluation of (Av,z)jim and _ of 
( Avz)um can thus be carried out graphically by inter- 
preting Eq. (48) in the way noted in Fig. 10b. The 
value of \; is then obtained directly from Eq. (47 

(C-1-1-3) Description of the velocities obtained afte 
the expansion at the leading edge takes place: Now con 
sider any point P," whatsoever, belonging to the char 
acteristic line o,') just determined and lying in the 
= R and also belonging to the characteristi: 
0 of the same plane 


plane Y 
line o,‘” passing through point 
it is considered to be close enough to the leading edge 
to make the are Py Py, of the line o;’ indistinguish 
able from its tangent at Py. This is the same situa 
tion already met with in Section (B-1-2), and the veloc 
ity at Po can be obtained by means of the construc 
tion therein described (Figs. lla and 11b). 

Since it is legitimate to set fy = 0 in the plane under 
consideration because of the assumption made at the 
beginning of Section (C-1-1), then knowing the velocity 
at Po‘ makes it possible to calculate \; by aid of Eq. 
(39), and therefore the tangent to a, at Po may be 
drawn in. 

In the plane Y = R + A, consider the point of inter 
section QV) that is formed when the tangent to o,"" at 
the leading edge, and already drawn in [see Sectiot 
(C-1-2)], meets the plane Z 
the point Po“. Let Qo be the point of o, 
on the same characteristic line 0,’ that passes through 
and it 


const. passing through 


which lies 


Qo". Then the velocity components at Qo"! 
the plane } = R + A may be determined by means 
of Eqs. (42’) and (43’). Meanwhile, the difference 
between the values of vy at the points Q,"? and at the 
leading edge A [or at the points Qo? and Qo] may be 
computed through use of the fourth of Eqs. (42), or the 
fourth of Eqs. (43); consequently, the value of zy at 
Qo? is determined. 

Let the equation of the surface }°” in the neighbor- 
Zp 


hood of the plane Z be written in the form 


A Zp A, } Tt VA Zi A 


then there exist the relationships 


and 


(fy)z ZPypD = ZP, pA (J 


Now, since A,°(Y) is known, substitution of this 
value in Eq. (39) allows the calculation of \, from tt 
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without ado, and then the tangent to o;? at Qo may 
be traced in. 

Thus the repeated application of this procedure pro- 
duces the entire surface =“ and the velocities at its 
various points. The determination of the >“ in the 
plane Y = R is then undertaken, but this does not re- 
quire anything other than the repetition of the method 
just described. 

(C-1-2) Determination of the Field Along the Top Sur- 
face of the Wing.—-In the plane Y = R, let D be the 
point at which the characteristic line o;’, passing 
through Py” of the line o, and also sufficiently close 
to the leading edge, intersects the plane of the wing 
(Fig. 12a). The change in velocity along the element 
P\™D satisfies Eq. (43’); on the other hand, the 
velocity at D has a known direction, and on this ac- 
count the graphical construction, shown in Fig. 12b, 
is evolved for determination of the velocity. It is 
thus possible in this way to draw in the tangent to the 
characteristic line o; passing through D by letting 
fy-er = 0 in Eq. (39) and by using the values just 
obtained for v, and gz. 

In the plane Y = R + A, let D’ be the point on the 
wing which lies in this plane and has the same abscissa 
as D, and let Qo™ be the point of the line o,;% which 
belongs to the same characteristic line o’ passing 
through D’. Since the velocity at Qo” is known (or 
can be determined immediately by interpolation be- 
tween the values of the velocities computed along 0," 
by means of the above-mentioned procedure) and 
since at the same time the direction at D’ of the veloc- 
ity component in the plane under consideration is 
known, then the magnitude of- this component can be 
determined by the usual graphical construction (Figs. 
l2e and 12d). Furthermore, because of the irrota- 
tionality of the motion, we know Ovy/Ox = 0v,/OY, 
and, therefore, if A denotes the leading edge of the 
wing, it follows that 


lhe value of (vy)p may be obtained easily from this 
relation. 
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graph for the field over the upper surface of the wing 


Eq. (39) gives the value of \; if one assumes for the 
velocity components the values just determined at D’ 
and also sets (iv)gue = (0 and (fz)se0 \, therein. 
The calculation of this quantity thus makes it possible 
to draw in the tangent to o; at D’. 

In this way one may obtain the velocities at all the 
points in the plane of the wing which belong to the 
plane x = Xp, and proceeding in the usual way one 
may also determine the entire characteristic surface > 
which intersects the wing along the straight line Z = 0, 
Xx = Xp, and one may calculate the velocities at the 
points of this surface. 

Therefore the complete determination of the flow 
over the upper surface requires but the repeated appli 
cation of the procedure just described. 


C-2, Description of the Field Existing on the Bottom Surface of the Wing 


C-2-1) Determination of the Shock Wave for the Flow Field in Contact with the Lower Surface. 


-The compression 


in the fluid stream which touches the lower surface is propagated by means of a shock wave that originates at the 


leading edge of the wing. 
pressed in the form: 


Let the equation of the shock wave, in the neighborhood of the leading edge, be ex- 


x = Zg(Y) 1) 


Since the direction cosines of the normal, n, to the wave surface are, for x = Z = 0: 


1 


Vi1+ g? 


cos (n, x) = 


cos (n, Y) = 0; 


cos (n, Z) = , 
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then the condition that the velocity be normal to the wave front leads us to write 


(0) _ a 9),,(0) 
er a 3 vy 


l 0 


’ denote the values of 


0) o..(0) o_(0 


where we again let v,“", vy“, vz 


denote the values of the same components after the shock. 


= Ur » — i: Vz = vz" + gil 


Pd (0). a 
7 — 4 >. = 
v 


a mies Cy , ez 


After the shock the stream must be tangent to the wing, and therefore it must be true that #z = 0; 
Consider, in addition, that the velocity components normal to the wave are 


H = —vz/g. 
[v, i/] + g?| oe 


before the shock and 


7 1, 


! 


ip eae a 


7 
the velocity components before the shock and 2,, dy, i, 
From this set of equalities we may write 


consequently, 


[v2(g)/W1 + g?] 


[v,, V1 + g*] — [/7/V/1 + g?] 


behind the shock, while the tangential component is expressed as 


vf = flor + for? + fez} sin? 


= [vy + 


[v, |?g? + [vz |? + 27, 7g gz 
1 + g? 


provided y denotes the angle between the normal, 7, to the shock and the direction of the upstream velocity. 


Then, by utilization of the equation of the shock, it follows that 


v,0 Vz, g v, vz 
| ; , | / , . | : 
Vicg Vi¢@PIlVic+g fVicg 


that 1s, 
k— 1 


vs — v2 glloee + v2) _ 
k+1 


gil — [vy ]}*} 


This relationship is equivalent to the equation defin- 
ing g for a two-dimensional stream that has, ahead of 
the shock, the components 

gS" V1 — [vy |? and vz" V1 
Remember now that 


vy, = V/V, and vz = Vz/V, 
and set 

i tte 5 tal A tad (51) 
so that 


0 


v,(° V, Uz Vz 


and 
V1 — [vy]? 


/ 70) Nos V* 


It is apparent then that, in any plane Y = const., the 
trace of the shock wave at the leading edge can be de- 
termined as though the upstream flow were two-dimen- 
sional and with components equal to the actual com- 
ponents V, and Vz but with the limiting velocity 
of the stream less than the real one by the amount indi- 
cated in Eq. (51). 

(C-2-2) Determination of the Field Along the Lower 
Surface.—(C-2-2-1) Construction of the series of charac- 


[oe }g? + loz P + 2ve 02g | 
1+ g? 


0 J 
2) 50 


J _ y,,(0) ]2 — 


v2 Pg? + [oz P + 20,7! 
11 — [oy P31 + g*) 


teristic surfaces that emanate from spanwise lines. aft of 
the leading edge: Inthe plane Y = R, let So“ be a point 
on the trace of the shock wave and situated so close to 
the leading edge that it is legitimate to consider S, 

as belonging to the tangent to this trace, which was 
described in Section (C-2-1). Now consider the 
characteristic line o,’, passing through S)‘"; its slope 
can be determined by means of Eq. (39) because 
fy = (0 and the velocity at Sy“ Let D, 
denote the intersection of o; with the plane of the wing 
(Fig. 13a). Then the velocity change along S)\"D, 
satisfies Eq. (43’), and since vz must also be zero at Dy, 


is known. 
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Fics. 13a-13b. Determination of velocity on lower surface of 
wing behind the compression shock. 


ere 


O8 ERI eo 


yr, Vz 


ntly, 


city. 


0) 


ift of 
yoint 
se to 
So 
was 

the 
slope 
“ause 
t D; 
wing 
OD, 
t D;, 


E y-R 


ce of 





INTERFERENCE 





THE HODOGRAPH FOR THE PLANE Y«R 


THE PLANE Y¥*R 





liv, 7 3,. 

fo) | | 

Jan 

V"l vevociry ar SO aHeaD OF 
THE 


”, vevocrry at S\leHino 
Tre: SHOCK 





V 








Determination of velocity on the shock front 
for lower surface of wing. 


l4a—14b. 


FIGS 


the velocity at Dy, is found with little difficulty by 
means of the construction given in Fig. 13b, for ex- 
ample. 

Now that the value of the velocity is known at D,, 
it is possible to compute the value of is = \, by means 
of Eq. (39), wherein fy = 0 is still to be employed; 
consequently, it is possible to draw in the tangent to 
o\ through the point Dy. 

In the plane } = R + A (and in the other planes, 
) = const.), the same procedure as described above is 
followed in order to obtain the velocity at the points 
D,’, which lie on the same straight line, x = Xpp in the 
plane of the wing. At these points the tangent to 
o;'". which passes through them is traced in for each 
one of the planes. Proceeding along the line o,""’, it is 
possible to calculate the values of the velocity in each 
plane }’ = const. by means of the repeated application 
of the procedure already fully described, until finally 
the intersection of the surface =“ with the shock wave 
is reached. 

(C-2-2-2) Determination of the velocities at the inter- 
sections of the characteristic surfaces which begin on the 
wing surface and meet the shock wave from the leading 
edge: Let S, denote the intersection of the line o,"”, 
lying in the plane Y = R, with the trace of the shock 
wave. The determination of the velocity at S,“ and 
of the tangent to the trace of the shock wave at S,,” 
comes about as a result of the following considera- 
Along the are S,_;S, of o1’ (Fig. 14) the 
Since the veloci- 


tions. 
velocity change satisfies Eq. (42’). 
ties are known at all points of o;° which correspond 
to the S,_,, lying in the various planes Y = const., 
this Eq. (42’) furnishes the means of computing the 
velocity component Az; in the direction of o;’ and 
through S,_,;. This component must be added, to- 


gether with its normal component, to the velocity at 


BETW 
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S,-1° in order to obtain the velocity at S,. On the 
other hand, at S,, downstream of the shock, the 
velocity change produced by the shock must satisfy the 
equation of the shock wave itself. If the vector defin- 
ing the velocity ahead of the shock, at S,,, is drawn in 
in the hodograph plane, together with the correspond- 
ing shock polar, then the construction shown in Fig. 
14b is readily available for determination of the velocity 
behind the shock at S,,. 

In the plane Y = R + A, let 7, be the intersec- 
tion of the line o,“” with the trace of the shock wave. 
The points S,, T,,, ete., which are obtained in the 
various planes Y = const., constitute the locus of a line. 
Starting at this line, the equation of the shock wave 


may be expressed in the form 
x= A’(Y)g(VY) + (Z — A’)g,( VY) (52) 
wherein A’( J’) is the coordinate relative to the Z-axis 
of the points S,, T,,, ete. The direction param- 
eters for the normal to the shock surface at the points 
indicated become: 
APP. 
1; —A’g ~ A¢-+ A’g, = —A’g; -g 
respectively, provided it is assumed that A’(g — g;) 
may be neglected on the basis that it is of second order 
in magnitude in comparison with A’g. 

In each plane Y = const., assume, in addition, that 
the x-axis is the straight line parallel to the velocity 
component at 7,“ ahead of the shock (this point of 
view is convenient in order to reduce the formal com- 
plication of the computation; because of this, the 
values of A’ and of g will be modified in an obvious 
way with respect to their corresponding values in the 
system of axes which has been used until now). Then, 
employing the same notation as adopted in Section 
(C-2-1) and with the understanding that the velocity 
components are now referred to the system of axes 
which was just described, it must be true that, as a 
result of the condition that the velocity change pro- 
duced by the shock has to be perpendicular to the 


wave, 
v (0) — v dy —_ vy Vz 
Zz r A 
] — A’g £1 


from which we obtain 


v, =7,% —H; dy = vy + A’gH; v2 = gH (53) 


If such terms as A’’g?H and [vy |?A’’g? are neglected, then use of the shock equation produces the rela- 


tion: 

7,0 vy A’p | Vz vy A’g 

; x —_ A J = — f a ren 
Vi¢g? Vit¢gilVi¢g? Vit ge 


Uz21 ; | ‘= 
Vi + gi” 
bm J f, 7 lop} : g,°[v, |? + 20, vy A’g) (54) 
k+1\ + gi pis 
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If the velocities are referred to the reduced limiting velocity defined by Eq. (51) instead of being referred to the 
actual limiting velocity, so that we therefore let 





vy, vz 
a a ee - = 9,*- t 
; ae ig = Ug , CU 
V1 — [vy]? vi = [vy }? 
then Eq. (54) becomes 
k— 1 k— 1 l 
no + ee a ae 2 *\2,.2 Joy Fay ¥ ta o 
02°02" — Uz, £02," = (1 + £1") — (v, "£1 + 20," 0y" 4 g 
" k+1 k+1 k+1 
wherein the term 0z*vy* A’? has now been neglected and we have let 
v,*vy* A’g = Vz *vy* A’ 
But from Eq. (53) we have 
£1 = 02/H = iz/|v, — 3,) = 0z*/(v,* — 0,*) 


so that, upon solving for #z**, it is found that 


vz** = (y,* — 9,*)? iF = = ." 5,” + 2v,*vy* am. @) | /Tee v,* - a aa _ien : { 55) 
Le +1 k+1.° 41 £4153 


Substitution of the value vy* = 0 (or of g = 0) into this equation reduces it to the well-known relationship that 
1 I 


defines the shock strophoid. 

Eq. (55) describes the new shock polar in the arbitrary general plane Y = const. Once this polar has been 
drawn in in the respective hodograph it permits the determination of the velocities at the points 7, and of 
the traces of the shock wave at these points, by means of a procedure similar to that discussed in the case of the 
hodograph for the plane Y = &. In this way the complete determination of the flow along the bottom surface 
requires only the repeated application of the procedure described in the preceding sections. 


(D) DETERMINATION OF THE FIELD FOR THE PORTION 
OF THE Bopy TuatT Is AFFECTED BY THE WING 


D-0) General Equations 

In order to study the flow near the fuselage in the region where the disturbances produced by the wing are felt, 
it is still convenient to refer the velocity potential of the field to a system of cylindrical coordinates, as was done 
in Section (B). Now solve the equation of the characteristic surfaces for x, and write the equation so obtained in 
the form x = f(y, 0, b) through use of such a cylindrical reference system, so that then 


: Pes Pe as a , a 
Gn = PAyafy oad 2a43 ( i oh 223 ( fit, 4 Asa ( - fy? 2 Aeofy” = () na) 
3 y y y* . 


The intersections of the set of characteristic surfaces > with the cylinders y = const. and with the planes @ = 
The fundamental equations that 


const. produce the characteristic lines that are still to be denoted by o; and oz. 
define the law of variation of the velocity along these characteristic lines now become 


. f = N; hy = py 


33 [d (yup) dv, dv, dv, 1 d(yuy Vy 

E ( de T Ha dé )| r E (201s Gash) 1 Os dy se vy dy | r y . 7 57 
d (yup) dv, dv, , dv, 
| dy ** dy | ~ G bil dé ). | 


At the same time, the intersections of the set of characteristic surfaces ¥’ with the cylinders y 
planes 6 = const. produce the characteristic lines denoted by o;’ and a’, and for these lines, the following hold 


ao (dome ie, ioe. | ie, 1\ d(yn ° 
E (‘ (y 0) + M1 : ) bh f (2aj2 = Ad) + 22 7 s + 223 () = am | + a = 0 
y? dé dO J |e,’ dy dy yi dy jay Y =x) 
| 2% 4 a | - (2 +> “*) | ; 
dy = dy}, \d0 ~ do/,, 
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INTERFERENCE 


BETWEEN 


WING AND BODY 


fhe quantities yw; and ye in Eqs. (57) and (58) are the two roots of Eq. (56) which result when f, is set equal to 


and the equation is solved for fg. 


Eqs. (57) and (58) are valid over the whole region that corresponds, through 


Eq. (56), to the case where the two values of /,, obtained for a given value X of f,, are real and distinct; this cer 


tainly happens in the neighborhood of the plane @ = 0 (or 6 = =). 
Analogously to what was done in Sections (B-0-2) and (C-0-2), the third of Eqs. (57) and (58) may be replaced 


‘ 1d; . v . , =F hd 
E ( | + | | sin (72, x) + (*) sin (72, x) = 0 se 
y \dé /I,, ve ¥ 
Fase (dv v SN’ 
(4)] +[] miner + amine 
y dé io? baat y 


by: 


Here 7; and 72 are, respectively, the tangents to o, and o,’ at the point P, in the neighborhood of which the 
In addition, dv; /ydé defines the change along o; of the component of veloc 


(57’) and (58’) are being considered. 


Eys. 


ity in the direction of the tangent to o;’; likewise, dv, /vd6 defines the change along o;’ of the component of veloc 


ity in the direction of the tangent to o,. 


D-1) Disturbed Flow over the Body Emanating from the 

Upper Side of the Wing 

The perturbations produced by the wing will first be 
felt starting at the characteristic surface =“; just as 
was done in Section (C-1-1-1), this surface is determined 
first of all. In the neighborhood of the leading edge 
of the wing and therefore in the neighborhood of 6 = 0, 
let the equation of ©‘° be expressed in the form 


x = Oho (y) (59) 


Then we have 


d= (feng = 0; wi = (feleno = ho(y) 
which, when substituted into Eq. (56), leads to the de- 
termination of hy(y). It is readily perceived that the 
quantity /o?(y) coincides with Ao“(}) as defined in 
Section (C-1-1-1). Now consider the various cylindri- 
cal surfaces y = const. (in particular, the one y = R), 
or, even better, let us give our attention to the planes 
into which these cylindrical surfaces can be developed. 
Then, in view of the above, the tangents to the lines 
7” can be drawn, where o; are the intersection of 
these planes with the surface S. 

The surface [Eq. (59)] is then cut by the plane 6 = 
4. In the neighborhood of the line of intersection x = 
Vi (y), 0 = A, it is assumed that the equation of 


< is given by: 


x = Aho (y) + (6 — ADA, (y) 


from which it now follows that 


A = (fie=a 


= Aho 
Mr = (folona = n(y) 

Thus Eq. (56) affords the means of obtaining /y(y) 
when the values just given for /, and fs are substituted 
init. Then it is possible to trace in all the o, by 


Proceeding step by step over the cylindrical surfaces 


f 


or, better, upon the planes into which the latter are de 
veloped. 

With respect to what was done in Section (C-1), 
therefore, no modification arises in the procedure used 
in application of Eqs. (57) and (58) and Eqs. (57’) and 
(58’) to the determination of the characteristic lines 
o, and o;’ and oe and oa’, as well as of the related hodo- 
graphs. The procedure is now carried out on the planes 
into which the cylindrical surfaces y = const. are de 
veloped in the same way as was done in Section (C) for 
the planes Y = const.; likewise, the planes 6 = const. 
now correspond to the planes Z = const. 

Because of the considerations evolved in Section (C) 
it is indeed sufficient to consider the flow only in the 
R. When 


one restricts himself therefore to so consider the motion, 


neighborhood of the cylindrical surface y = 


the equations can be simplified, in actual calculation, 
because v, = O for y = R; so we may neglect the term 
v,/y, and may let da = | and a2; = a3 = O therein. 
One must then consider a set of cylindrical surfaces, 
y = const., which are close to one another in order that 
such simplifications would be legitimate for them. 
Moreover, they must be of sufficient number to allow 
the determination of the flow over the entire portion 


of the body which is of interest. 


(D-2) Disturbed Flow over the Body Emanating from the 

Under Side of the Wing 

Observe that on the cylinder y = R, which is a stream 
surface, the shock wave is normal thereto because the 
shock must not produce any change of the velocity 
component v,. Therefore, in the neighborhood of this 
surface the equation of the shock front may be ex- 
pressed as 


9 


x = go(O) + gi(O)[(y — R)?] (60) 


The following conditions serve to determine go(@) and 
gi(0): The shock wave has an intersection (trace) 


on the cylinder y = const. The tangent to this trace 
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at the plane @ = 0 is found easily and directly by means 
of the deductions made in Section (C-2). ° 

Consequently, one proceeds, in the manner that was 
there described, by working in the planes into which 
these cylinders have been developed; and through 
utilization of Eqs. (57) and (58) and of Eqs. (57’) and 
(58’), the velocity at the points D,; is calculated. The 
points D, belong to the plane of the wing and, of course, 
are close to the leading edge. 

Now, as long as the characteristics o, that originate 
at these points D, do not intersect the trace of the 
shock wave on the respective cylinders (that is, on the 
related developed planes), the velocity component per- 
taining to each of these cylinders has the direction of 
the x-axis, downstream of the shock front. Then, 
because the velocity upstream of the shock is known, 
the change in direction that the shock must produce is 
consequently also known. 

It is possible, therefore, to determine the trace of the 
shock which is capable of producing the required 
change through evaluation of the shock polar (which 
coincides with the shock polar based on the local veloc- 
ity, because v, is zero or negligible, in the region of the 
field under consideration). The characteristics o, 
through the points D, of the wing contour, intersect 
the traces of the shock which were just described. 
Again, through employing the method developed in 
(C-2-2-2), the velocities at these points of intersection 
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are calculated (for the case where the cylinders are de- 
veloped into planes). One may assume that the shock 


polar is the strophoid that corresponds to the local 


velocity, in view of the fact that v, is negligible; thus 
the above-mentioned method is simplified. 
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Erratum 


The title of the article by Adolf Busemann in the June, 1949, issue of the JOURNAL (Vol. 16, No 


The correct title is: 


“The Problem of Drag at High Subsonic Speeds”’ 


or Busemann, A., Driicke auf 
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Performance Analysis of Centrifugal 
Compressors 


* SHAO-PAN LIANG* 
Chekiang University, Hangchow, China 


ABSTRACT 


Elementary analysis of centrifugal compressors indicates that 
the energy received by each unit weight of gas compressed is 
merely a function of impeller tip speed. This, as shown by experi- 
mental results, is far from true. The pressure of delivery falls 
off with the increase of quantity of discharge, and the efficiency 
of compression is also affected. It is the purpose of this paper to 
point out that, based on force and energy calculations, perform 


ince can be predicted with reasonable accuracy. 


INTRODUCTION 


pe DEVELOPMENTS OF POWER GAS TURBINES and 
jet engines have opened a new field of application 
jor compressors of the centrifugal type. Following 
Moss's analysis, the energy input per unit weight of gas 


1S 


W = SU?/g (1) 


where U is the tip speed of the impeller; g, the gravita- 
tional acceleration; and S, the ratio of tangential veloc- 
ity of the gas leaving the impeller to the impeller tip 
speed. The discharge pressure depends on S and yn and 
the efficiency of compression and has not been theoreti- 
In this analysis, the author has fol- 
momentum and 


cally surveyed. 
the 
energy considerations. 


conventional method of 


Such relations alone, however, 


lowed 
do not yield fruitful results. The friction loss in the 
impeller is an unknown factor, which makes the process 
deviate from the isentropic one. By treating the proc- 
ess as a combination of losing energy and adding heat, 


an extra equation is obtained. 


ENERGY RECEIVED BY AND FORCES ACTING ON AN 
ELEMENTARY MASS OF AIR 


Consider an elementary mass of air between two 
adjacent blades and at radius r on an impeller (Fig. 1). 
The angular momentum about the axis of the impeller 
is Sr°wpA dr/g, where w is the angular velocity of the 
impeller, p is the density of the air in question, and A 
is the cross-sectional area of airflow. Since the angular 
momentum of the air is zero before entering the impeller, 
this change of momentum is associated with an increase 
of energy, which can be easily seen to be Sr*w?/g per 
unit weight of air. This energy may take either the 
form of sensible heat or stored as mechanical energy, 
and the total sum is 
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(Srw)? wv ee 
a + +t C,(7 


22 ~g 


T (2) 
in which v is the relative radial velocity of air in the 
impeller; C, is the specific heat of air at constant pres- 
sure, expressed in mechanical units; 7 is the absolute 
temperature of air in the impeller; and 7% is the abso 
lute temperature of atmospheric air. 

The forces acting on the elementary mass of air are 
(a) the pressure of the neighboring air and of the im 
peller blades, and (b) the friction forces exerted by the 
impeller. 

Taking the above-mentioned air as a free body, the 
summation of forces along radial direction is 


pu* 


Ad 
— dp o_o ¥2 p ar i » 
-£ 


bdr = 
g 

where f is the coefficient of friction of the impeller, / is 
the length of circumference of the area concerned, and 
a, is the acceleration of air in radial direction. 

As the motion of the air is composed of a circular 
motion with tangential velocity Srw and a radial motion 
with velocity v, the radial component of acceleration is 


a, = v(dv/dr) — S?rw? (4) 
Hence, 
v Ad ly 
—Adp — ge ldr = P ~4 (0 - rs) (5) 
” 2Qg gR1 dr 


in which R is the gas constant of air. 

Since the frictional force is, in general, extremely 
small compared with the pressure gradient, the equa- 
tion can be further simplified by omitting the friction 
term, with the result 


R dP ] (: dv 


 pdR~ eT 


— Sry? (6) 
dr 


Eliminating v between Eqs. (2) and (6), we have 


R dp  --- d 
p dr — 20T ldr 


_ 2Sr°w* 


— §*r2w* — 


2C,(T — To) | — 2sratt (7) 
By assuming that S is independent of 7, the equation 
may be simplified to 
R dp l 
p dr a gT 


7 _ aT 
2S(1 — S)rw? — gC, 
dr 
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Fic. 1. Forces acting on air at radius r. 


This gives the relation between pressure and tempera- 
ture at different radii, and another equation is required 
to obtain a solution. 


THE TRANSFORMATION OF ENERGY DUE TO FRICTION 


In order to determine the pressure variation, we must 
consider the effect of friction loss. With the assumption 
of no heat transfer from outside, any loss of energy due 
to friction will again be absorbed by the gas. There- 
fore the process is practically a transformation of me- 
chanical energy into heat energy, with the total amount 
unchanged. 

Since the cross-sectional area of airflow is roughly of 
a rectangular shape, the length of circumference may be 
expressed as 


1 = 2[(2ar/n) + (nA/2zr)| (9) 


in which 7 is the total number of blades. 

The energy loss per unit weight of gas is directly pro- 
portional to the velocity head (relative velocity) and 
lateral area, while inversely proportional to the cross- 
sectional area of flow, so the loss due to impeller fric- 


tion may be expressed as 


— §)*r?w? 2ar nA\ . 
2 + ar 
22 n 2nr 


The friction due to the compressor casing, being less in 
magnitude for normal rate of flow, may be neglected. 
Thus, the entropy increase due to this energy trans- 


vv? + (1 
f 


formation is, from gas laws, 


2fv? + (1 — S)*r?w? 


7 p 1 2¢ 
2arr n 
( + ) dr (10) 
nA 2nrr 


This relationship between friction and change of 
state helps in determining the effectiveness of com- 
pression. As can be readily seen, the less the magnitude 
of f, the nearer will be the process to the isentropic 


i 


one. 
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Eliminating v between Eqs. (2) and (10), we have 


CG,dF Rdp 2 & a | 
ee - = = ( ) 7 —_ To) 
T dr pd TL 2% " " 


(2x 4 n ) ' 
| 
. nA 2nrr 


Substitution of Eq. (8) into Eq. (11) results in 


S(1 — S)rw? 


rw? 
: =f - 
g 2¢ 


GAY — ro | x 


2rr n 
+ { 12) 
nA Qnr 


which gives the absolute temperature at different radii. 
With this known, the pressure of gas can be found from 
Eq. (8), while the radial velocity v is found from the 
combination of Eqs. (4) and (12)—-namely, 


— S)rw? 


v? + (1 — s?)r?w? 
g <£ 
(2= _* ) ” 
nA 2nr P 
After the computation of pressure and temperature, the 


quantity of discharge becomes extremely simple as 
indicated by the equation of continuity 


S(O 


gq = punA 14) 


where gq is the weight of discharge per unit time. 


SIMPLIFICATION OF RESULTS 


The above equations, while based on _ theoretical 
analysis, are rather complicated in actual computations. 
They can be further simplified if the number of im- 
peller blades is small. In such case, the term /2zr is 
sinall in comparison to the term 2x7/nA and may be 
dropped out. Thus we have the approximate equa 


tions 
S(1 — S)w? | r2w? ons _. oe ; 
=f — C,(T — To) (15 
g 2g nA 
and 
S(1 — S)w? = (fax/nA)[v? + (1 — S)?r?w?] (16 


instead of Eqs. (12) and (13). 
Since A is practically constant, we have dA/dr 
0, so that differentiation of Eq. (15) gives 


CAleT ar) = ra?/g (li 
Combining with Eq. (8) gives 
R(dp/p) = ZC,(dT/T) (18 


in which Z? = S? + (1 — S)*. This holds good for the 
entire path through the impeller. 

For the complete compressor, detailed information 
regarding the entry and diffuser section is required to 
obtain a definite solution. By assuming a constant 
temperature-pressure relationship throughout the unit. 
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the problem can be further simplified. Thus, the equa- 


tion 
ps/po = (T/T) ~» (19) 
gives the pressure ratio and the expression 
n = [(T3/To)’ — 1]/((13/To) — 1] (20) 
gives the adiabatic efficiency of compression. (The 


symbol y denotes the ratio of specific heats of air.) 
The final temperature 7; may be determined from the 
law of conservation of energy as 


T; = Ty) + (SU?/gC,) (21) 


From the above statements, the density ratio at 


radius 7; is 
Dp) j : oe l ("*) 
-= 41+ x 
Po ' 2 Co 
SQL — S)\1 — hk) N77? 
E _ \ - 3 _ aa (22) 


in which p; is the density of air at radius 7, py is the 
atmospheric density, 7; is the radius of impeller eye, 
his the ratio of diameters of impeller hub and impeller 
eye, and ¢ is the sonic velocity in the atmospheric air. 
The rate of flow is 
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q = pwrr,?(1 —_ h?) x 


. . 9 I/s 
S(1 — S)(1 — &*) ee 
ee — (1 — S)? (23) 
J 
The specific volume-rate-of-flow, as commonly used in 
American practice is 
Q T pi Nw 


D/6 Vm D? 
S(1 — S)(1 — h?) _ 
: _ - — =~ {h 5)| (24) 


in which Q is the total volume of airflow per unit time, 
D is the impeller tip diameter, and @ is the ratio of 7 
(519°F. is 


(1 — h?) X 


to the assigned reference temperature. 
commonly used.) 


EMPIRICAL LIMIT OF APPLICATION 


From the above equations, it is evident that S is 
determined for given values of Q and f, at any speed. 
The pressure ratio, as shown in Eq. (19), is also fixed 
under such conditions. With the value of f properly 
chosen, theoretical results will check fairly well with the 
experimental ones. In Fig. 2, the data of Campbell and 
Talbert? are taken for comparison. The equation 
loses its validity when maximum pressure is attained, 
and this happens in the neighborhood of S = 0.9 for all 
speeds. Further reduction of quantity would theoreti- 
cally give still higher pressures, while experimental 
data show otherwise. The explanation for such dis- 
crepancy is that friction due to impeller casing was 
neglected in the above derivations. When the quantity 
of flow is small, the heat generated becomes greater and 
greater for each unit weight of air. Consequently, the 
temperature rises higher and higher and efficiency is 
lowered. The quick fall of efficiency in this range leads 
to the phenomenon of surge, as may be seen in prac- 
tice. By assuming that S is constant after reaching 
0.9, sufficiently close results can be obtained up to the 
surge point. Since operation in this range is not stable, 
only the right portions of the curves are of practical 
interest. 
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and discussions of papers published in the JOURNAL will be 
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Transient Phenomena at Sonic Speed 


N. Rott 

Institute for Aerodynamics, Swiss Federal Institute of Technology, 
Ziirich 

March 21, 1949 


B: APPLICATION of the linearized theory to the calculation of 
transient loads of an airfoil of infinite span and of zero 
thickness, finite results are obtained even at sonic speed, as in the 
case of harmonic oscillations.! With the latter solution as a 
starting point, the same methods may be applied as in the super- 
sonic case.2 (Further references are given in reference 2.) As 
an example, the response of an airfoil to a sharp-edged gust will 
be calculated for M = 1. 
If the velocity distribution along the airfoil in the plane y = 0, 
leading edge at x = 0, is given by a wave proceeding with the 
sonic velocity a, 


ikcat ) 
vo = —vze' . " (1) 
the corresponding pressure (on the upper side of the airfoil) at 
sonic speed may be calculated from Eq. (5) of reference 1: 
etka (r 2)) 


Po = —par(2riKx) (2) 


fo represent boundary conditions for an airfoil in a sharp gust 


with constant upward velocity 7, the step-function is obtained 
from Eq. (1) by the following operation :* 


No proof will be 























| [4] 
Tt 
j= | 
ic |7 
Fic. 1. OCz,/Oag vs. Ut/C. 
pe = (pav/r)V (2at/x) — 1, x < 2at (5a) 
= (, x > 2at (5b) 


The same result appears also as the limiting case of the super- 


sonic solution.?, According to Eq. (16) of reference 2, the pres- 


sure pg in a stream of velocity U = Ma (M > 1) is 
, M : 
pe = —pav / pe * x < (M — Il)at (6a) 
M y¥ 
= —pav > »(M — lhat<x< (M+ ljat_ (6b) 
V M*-—1¢f 
= 0, x > (M+ Dat (6c) 
where 
cosy = M — (M? — 1l)at/x (6d) 
For M — 1, cos Y — 1; then ¥ must be small so that cos y = 


1 — (2/2) or from Eq. (6d) 
y= VM — [CM + D(at/x) — YJ 


Introducing this limiting value in Eq. (6b), Eq. (5a) follows for 
M = 1, 

For the lift coefficient of an airfoil, with chord c, due to a gust 
at sonic speed, the following expressions are obtained from Eq 


(5) (introducing ag = v/a): 
l “FO wadK #4 
'%G= lim : -=-v7, x<al (3a) Oc, /O0ag = 4at/c, at < c/2 (7a) 
2r B 0 — B + ik : 
4at o + sin o * = 
=0, x<at (3b) - » arcfZ (7b) 
tf us 
The corresponding pressure is given by ee 
p= 1 — (c/at (7 
. - ° > be . => _ *) 
be = (1/2r) lim pod K/(8 + 1K) eicid ede ” 
es =a For large values of at/c, 
or introducing po from Eq. (2) a — ee 
' & polr 1 0c, /Oag (8 n)V 2at/c (7d) 
pav 1 t e* (2-2/2) . 
bg = - 7 lim 7 dK (A) In Fig. 1, 0c9,/agq is plotted, together with the corresponding 
Gane 9 ‘Kk Kk . . ‘ ° P -_ 
V 2nx <4 g = 0 2 (8 + iK)Vik supersonic values given by Miles? and the solution of Kiissner 
If and von Karman-Sears for M = 0. 
ea The distance d of the center of pressure from the leading edge 
/ . ° 
1(s) = (1/2) lim J e’’* dK/(B + iK)V ik is given by 
eee d/c = at/2c, at < c/2 (8a) 
then the inverse (or Carson’s) integral equation for A (s) is? : 
re o-— Sing cos¢ 
x = (at/2c) re a at > c/2 (8b) 
ss * o sin 
/VK=K f e~"*A(s) ds ’ 
0 For at/c > ~,d/c > 1/3. 
Which is satisfied by Sudden changes of the angle of attack may be investigated 
o/ Jf at sonic speed by the same methods. If the disturbance remains 
» = )/ . e . . . . 
A(s) = (2/V x)Vs stationary after the sudden change, it is characteristic for all 
Applied to Eq. (4), the pressure becomes solutions that the pressure grows finally witht. Naturally, 
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the validity of the results of linearized theory is at least doubtful 
for M = 1, but these finite solutions may suggest again the im- 


portance of nonstationary phenomena. 
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Oscillating Airfoils at 


Quasi-Stationary Thin Airfoil Theory 


John W. Miles 


Department of Engineering, University of California at Los Angeles 
March 17, 1949 


HE AERODYNAMIC DERIVATIVES for use in dynamic stability 

calculations are frequently calculated on the basis of steady 
flow, the argument being advanced that the motion is small so 
that nonstationary aerodynainic effects are negligible. This 
line of reasoning overlooks the fact that the damping forces are 
first order in frequency, and the quantities neglected may there 
fore be of the same order of magnitude as these forces. In par 
ticular, if @ is the velocity potential, the pressure is given by 
Newton’s law as 

Do Od Op 


_ —pU — —p— (1 
Pmt Pp 


b = j 
f Ox - ot 


The first term in this equation is the only one entering stationary 
airfoil theory, but the second is evidently first order in frequency. 

The correct lift and moment coefficients for quasi-stationary 
(i.e., correct to first order terms in frequency) motion of a two 
dimensional thin airfoil may be obtained as limiting cases of the 


flutter results.! In the case of incompressible flow the results are 


j rk 2 
: 1 — + ik| y — In xX 
{ k 
w(—cos ¢) 


(1 — cos vy) + ik sin® ¢ t dg + O(Rk*%ln k) (2 


[ (cos g=- cos 2 ¢) —ik(1 — cos ¢) X 


sin? ¢][w(—cos ¢)/U]de + 0 (k%ln k) (3) 


k = we 2U (4) 


where 7 is Euler’s constant (= 0.577216), w is the angular fre 


quency, c is the chord, U is the free stream velocity, ¢ = O(7r) 
specifies the leading (trailing) edge of the airfoil, and w is the 


local down-wash. If the displacement at a point x = —(c/2 
cos ¢ is 2(x)exp(twt), positive down, 
w(2x/c) = U(0z/O0x) + ws (5) 


If k is set equal to zero in Eqs. (2) and (3), they reduce to the 
For the case of a simple plunging 
: 1), w(x) is of order k, so that 


well-known Munk formulas.” 
motion of the airfoil (i.e., 2(x) 
the Munk formulas are correct to order k? In k, but in the case of 

is = [x + (c/4)]@}, Eqs. (2) 
and (3) yield 


Oc; j wk 3 2 1) ; 
= 2r4{1—- + ik + +7 — In | >+O(k7ln k 
da ] 2 2 k/ J 


OCm/Oa = —(wx/2)ik (4 


pitching about the quarter chord 


(6) 


whereas the results given by the Munk formula are 
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OC; Oa = 2x7(1 + ik) 6’) 
OCm/Oa = —(r/4)tk 7’ 


Since k will generally be quite small (<0.1), the failure of Eq. 


(6’) to include the logarithmic term is a matter of some impor- : 


tance, for the sign of the imaginary part of (0c:/0qa) is reversed 
While most of the damping in rotary motions of conventional 
aircraft arises from the plunging motion of the aerodynamic 
center, an appreciable part is due to rotation about this center— 
i.e., to the imaginary part of (0c;/Oa). If the distance from the 
airplane center of gravity to tail aerodynamic center is five times 
the mean tail chord, the decrease in total damping calculated on 
the basis of Eqs. (6) and (7) over that predicted by Eqs. (6’) 
and (7’) is approximately 38 per cent for k = 0.01. For smaller 
c.g.-tail chord ratios the decrease is, of course, greater. 
larly, the damping in pitch of a swept wing is reduced. The ef 
fect of this error in calculating stability derivatives would appear 
to be most severe in the case of fin damping of small, high-speed 
Thus, most modern jet aircraft tend to ‘‘snake,’”’ prob- 
(This 


Simi 


aircraft. 
ably due to reduced (damping) effectiveness of the fin. 





problem was evidently not so serious for lower speed aircraft for 
a number of reasons.) 

In conclusion, it should be remarked that compressibility 
effects in nonstationary flow are not correctly accounted for by 
simply introducing the Prandtl-Glauert factor (1 — M7?) 
since compressibility introduces additional phase lags. The 
writer has obtained a closed form solution for the Possio integral 
equation,* correct to terms of order k*ln k, and hopes to publish 
the results soon. The quasi-stationary problem for linearized, 
supersonic flow has been discussed previously.‘ 
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Maneuvering Horizontal Tail Loads 


Innes Bouton 
McDonnell Aircraft Corporation, St. Louis, Mo. 
March 30, 1949 


A SEMIRATIONAL EXPRESSION has been derived for use in ck 
termining the maneuvering horizontal tail loads for an 
airplane. This expression has the essential simplicity of the 
empirical methods now in use but includes in a rational fashion 
ill the factors that vary between different designs. It is derived 
from consideration of rational methods involving the airplan¢ 
motion resulting from a given elevator deflection variation with 
time but eliminates the complicated calculations based on an 
elevator deflection function that in itself is empirical. 
The net maneuvering tail load, P;, is defined as follows: 





P, = Cy. se Sw q ¢ a 61, 
lt lt 

where Cy... is the airplane pitching moment coefficient about 
the c.g. at the load factor in question, SB is the wing area, g 1s 
the dynamic pressure, Z is the wing mean aerodynamic chord, 
l, is the length from c.g. to tail center of pressure, 6 is the maxi 
mum or minimum angular acceleration of the airplane in pitch at 
the load factor associated with Cy... used, and Jy is the airplane t 
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moment of inertia in pitch. The first term in the expression is 
what is often called the balancing tail load, while the second term 
is the incremental tail load needed to produce the angular ac 
celeration, 6. 

The evaluation of the first term is straightforward, although it 
is often difficult to determine Cy, accurately. It is in the de- 
termination of the second term that current methods have a 
fundamental disagreement. 

Methods that specify an arbitrary value of 6 become irrational 
when statistical data, based on past configurations that were 
ll geometrically and aerodynamically similar, are applied to 
radically different configurations. For instance, there is a gen- 
eral trend in high-performance airplanes toward sweptback tails 
with relatively short chord elevators, which results in lower ele- 
yator effectiveness. Also the tail area and arm to the c.g. have 
increased less rapidly than the moment of inertia. Both of 
these trends result in a decrease in the ability of the horizontal 
tail to produce airplane angular acceleration. 

Methods that determine airplane motion from elevator motion 
satisfactorily account for the geometric and aerodynamic dif- 
ferences between airplanes, but the methods are tremendously 
complicated to use, depend on an arbitrary elevator motion, and 
are extremely critical to minor changes in the pattern of motion. 
Itis satisfactory at low speeds to assume that the elevator is moved 
as rapidly as possible to the maximum possible deflection as limited 
by elevator stops or stick force, held at that deflection as needed 
to develop the maximum load factor, and then reversed to some 
specified deflection as rapidly as possible. However, at high 
speeds, what is physically possible becomes psychologically im 
probable For instance, in a recent study on a high-speed air 
plane, extremely high angular accelerations were possible if the 
maximum load factor were developed in less than 0.2 sec. How 
ever, the stick motion is so critical that the continuance of up 
elevator motion for 0.02 sec. or three-quarters of one degree more 
than the critical value would result in exceeding the maximum 
load factor by 40 per cent. It is certain that the probability of 
such elevator manipulation is low. Pilots would be restrained in 
such maneuvers by the danger of exceeding the airplane design 
load factor 
current regulations by specifying an angular acceleration in 
speed or by specifying at high 


This reduced probability is recognized in some of the 


versely proportional to the air 
speeds an angular acceleration that is a given percentage of the 
low-speed value 

A formula for 6 which recognizes the variations between air 
planes but which uses an empirical constant to be determined 
statistically is presented below: 

: d5,\ (dCi/da), [da (> =) l 
4=B8B _ n - (2 
d; (dC; da), dé : fe So 4 
where @ is in radians per second per second and B is either +2.5 
or —3.5; dé,/dt is the maximum rate of elevator deflection in 
degrees per second (discussed further below); dCz./dqa is the lift 
curve slope of the tail and wing, respectively; da/dé is the ele 
vator “‘effectiveness’’; .S, is the horizontal tail area in square feet, 
!, in feet, and J, in slug feet squared as defined in Eq. (1); 7 is 
the maximum load factor reached during the maneuver; W/5S, 
is the airplane wing loading in pounds per square foot, and V, 
the airplane equivalent velocity in miles per hour during the 
maneuver. All the aerodynamic parameters are determined 
the Mach Number of the maneuver. 

Che +2.5 and —3.5 values for B are empirical and subject to 
revision, although they appear to be approximately correct for 
data presently available. The plus 2.5 value occurs at the be 
ginning of a maneuver, and the @ so determined should be used in 
Eq. (1) in combination with the Cy... associated with a load fac 
tor of 20 to 50 per cent of the maximum load factor reached in 
the maneuver. In ali cases, the load factor used in Eq. (2) in 
determining 6 is the maximum load factor encountered in the 


| Maneuver, even though the balancing part of the tail load is de 
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termined at some lesser load factor. The —3.5 value for B occurs 
at the time of maximum load factor, and the tail load angular ac- 
celeration increment so determined is combined with balancing 
tail load at the maximum load factor of the maneuver. 

In comparing experimental data with this formula, the maxi- 
mum rate of elevator deflection d5,/dt, which occurs during the 
period immediately preceding the peak value of 6, should be used 
in the expression. For computing a value of 6 for use in design, 
it is suggested that a value of d5,/dt of 75° per sec. be used, even 
The higher rates 
However, if 


though higher rates are physically possible. 
do not seem to result in appreciably higher loads. 
the elevator rate is restricted to some lesser value, as may be the 
case in power-operated controls, the lesser value should be used 
It is also possible that the 75° per sec. rate should be lower for 
extremely large airplanes, such as transports or bombers, but it is 
believed that the lower load factor to which these large airplanes 
are designed will reduce the @ in the proper proportion through 
the n term in the Eq. (2). 

To date, all experimental data compared with this formula 
have been within the limits of the +2.5 and —3.5 values for B, 
as have all calculated values made from elevator deflections con- 
sistent with the philosophy of this method. The writer would be 
interested in learning of values of loads or angular accelerations 
higher than those determined by the equations above 
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A Note on Detached Shock Waves 


E. V. Laitone 
University of California, Berkeley, Calif 
March 15, 1949 


I REFERENCE 1, a Taylor series expansion is obtained for the 
distance of a detached shock wave in front of a body of 
It would be of interest to compare these results 





revolution. 
with those of reference 2, wherein a simple approximation locates 
the distance (.S) of the normal portion of the detached shock wave 
ahead of the nose of the body (see Fig. 1). This approximation 
assumes that the subsonic velocity on the streamline along the 
axis of symmetry approaching the stagnation point downstream 
of the detached shock wave is the same as for incompressible flow 


with the proper boundary conditions. Because of symmetry, 
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the portion of the detached shock wave through which this 
streamline passed was a normal shock wave. For each free- 
stream supersonic Mach Number (M ,, ), this normal shock wave 
produces a known decrease to a subsonic velocity. This locates 
the normal portion of the detached shock wave as shown in Fig. 
1. Since the distance L of the sonic line could not be easily de 
termined, its effect was generally neglected in reference 2. 
Applying the method of reference 2 to a sharp nose slender 
body one obtains along the streamline leading to the stagnation 


point 
1 2 : u 1 L y'(x)y(x)dx 
. +(y-1)]/= =l1- — 
y+ 1] M3, u., 2), (S + x)? 
for a body of revolution and 
4 ; u l 
- (M, — 1) = =~ 1— r. BSE 
se u rV M,? — 1 
L y'(x)dx i 
. (S + x) : 
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for a two-dimensional body. These subsonic velocity variations 
include only the first order compressibility effects as obtained 
from linearized theory (see, for example, reference 3). An addi- 
tional approximation in Eq. (2) is that (17 — 1) is smallso that 
behind the normal shock wave 

Vi- 
as shown in reference 2. 

Now, in reference 4, it is shown that Eq. (1) can be rigorously 
applied to the limiting case of a very slender cone as the super 
sonic free-stream Mach Number approaches unity. For a con 
the sonic line is always formed near the shoulder of the con 
In the limiting case of a slender cone near a free-stream Mach 
Number of unity, this sonic line approaches a limiting valu 
given by a straight line normal to the cone axis and extending an 
unlimited distance from the shoulder of the cone (see Fig. 2 
Similarly, the limiting value approached by the detached shock 
wave is another straight line normal to the cone axis but an un 
limited distance ahead of the nose of the cone as indicated in 
Fig. 2. For this case the subsonic flow region is effectively of 
unlimited extent, and it is shown in reference 4 that the boundary 
conditions are satisfied by Eq. (1) for y(x) = @x and y’(x) = 
0, so that 
1 + (S/L) l 


) log SL _ 


6? 


e 1 


i+ (S/L)E 441” 


Say are 
9 
2 4 
-+(y-D]=1- (M 
M ,? y¥+1 
The additional approximation in Eq. (3) assumes that (WV. — 
is small so only its first power is retained (see reference 2). 
Then, for air, we have y = 1.400, so Eq. (3) may be written 
f 1 + (S/L) 1 ) 
M = 1 + 0.30 < log : - . 
) S/L 1 + (S/L) 
This variation of the detached shock wave location (s) with the 
free-stream Mach Number (./q is shown in Fig. 3 for various 
values of the semivertex angle (@) of a conical nose of length L 
The minimum values of the free-stream Mach Number for which 
i conical shock wave becomes attached to the conical nose art 


also shown in Fig. 3 by means of a small arrow. 
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Errata—Aerodynamic Forces on an Oscillating 
Flap at Supersonic Speeds 


John W. Miles 
Department of Engineering, University of California at Los Angeles 


April 13, 1949 


— FOLLOWING ERRORS have been found in my paper (JOUR 
NAL OF THE AERONAUTICAL SCIENCES, Vol. 15, No. %, pp 
565-568, September, 1948). 
a) The ordinates of Figs. 4 and 5 should be labelled [2/ (| 
c)]*k*Rsg and [2/(1 — c)]*k*/ng, respectively. 


b) Eq. (20) should read 
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( Eq. (22) should read 


Aaa k,M) + cA a(k, WV) 


1,a(k, VW = 


¢ 7 1+ < 
V(t.) 4 
2 
k, ) + 


= eat) 


~ 
_ 


Che numerical results given in the paper are rather inaccurate 
particularly those for R,a and J,«), and more accurate results are 
referenced in a recent report by Biot, et al.! 
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A Method for Accelerating the Convergence of 
an Iteration Procedure 


jabriel Isakson 

Aero-Elastic and Structures Research Laboratory, Massachusetts 
Institute of Technology, Cambridge, Mass 

April 4, 1949 


I MOST PROBLEMS in which an iteration or successive approxi 
mation procedure is used in order to obtain a solution, the 
rate of convergence decreases as the iteration proceeds. Conse 
quently, high orders of accuracy become increasingly difficult 
to achieve, and a large number of iterations may be required. 
In such problems a great deal of time and labor may be saved if 
the values in a small number of successive approximations can 
be used to predict the final correct value 

In the present method the assumption is made that the suc 
cessive values obtained in an iteration procedure follow closely 
the form of a simple exponential approaching its asymptote. 
hat is, if one value is considered as being the first in a series, the 
value obtained after the mth subsequent iteration is given by, 


x =X + Ae (1) 


vhere Y is the asymptotic value and A and a are constants. If 


ve, %;, and x. are three values obtained in successive iterations, we 


may write, 
y= ¥ +A 2 
vy, = XV + Ae (3) 
vo = XN + Ae (4 


Subtracting Eq. (4) from Eq. (3) gives, 
X} ~~ X= Ae~7(1 — « (@)) 
Dividing Eq. (6) by Eq. (5), 
(X, — X2)/(Xo — X11) = « ( 


Substituting Eq. (7) into Eq. (5) we obtain, 
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Xo — X 
Az - () 
1 — [(x, — x2)/(xo — &»)] 


Substituting Eq. (8) into Eq. (2) yields, finally, 


A = Xo - 


[ (x9 — x))2/(x9 — 2x, + X2)] (9) 


It is thus seen that the asymptotic value may be simply deter 
mined from three successive approximations. This asymptotic 
value will not, in general, be the true value being sought, but it 
will frequently lie much closer to it than the three values from 
which it was derived. 
obtained in this way may be used as the starting point for sub 


If greater accuracy is required, the value 


sequent iterations and further acceleration. 

It should be noted that the calculation involves small differ 
ences between successive approximate values, so that a fairly 
large number of significant figures must be carried if accuracy is 
to be achieved. 

The method has been applied successfully in the matrix itera 
tion procedure for determining the natural mode shapes and fre 
quencies of vibration of an elastic system. The convergence of 
the mode coefficients to their correct value is accelerated in this 
manner, and it has been found that in many cases two or more 
significant figures are gained in each acceleration cycle. It is 
frequently necessary to defer the acceleration till five or six itera 
tions have been performed, since the variation in the first two or 
three iterations is sometimes erratic. 

In some iteration solutions successive values oscillate about an 
asymptote. In such cases the acceleration cannot be applied 
directly to successive values but may be applied to alternate 
values. <A simple rough test for the applicability of the method 
is that the difference between the second and third values should 
be of the same sign and somewhat smaller than the difference be 
tween the first and second values in the sequence. 

It is interesting to note that Aiken’s 6? process,' developed spe 
cifically for application to matrix iteration in the solution of chat 


icteristic value problems, yields the same expression as Eq. (9) 
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On the Compressibility Correction Factor for 
Axially Symmetric Bodies 


W. H. Dorrance 
Research Associate in Aerodynamics, Aeronautical 
Center, University of Michigan, Ann Arbor, Mich. 


March 31, 1949 


Research 


—— PAPERS have appeared dealing with the problem of the 
application of the Prandtl-Glauert method to subsonic com 
pressible flow over three-dimensional bodies.'~*4 Attempts have 
been made using this notion to establish a compressibility cor 
rection factor for three-dimensional bodies similar to the familiar 
compressibility correction factor (1 — M?%)-'/? applied to two 
dimensional bodies. It is generally agreed that the Prandtl 
Glauert method will not yield a universal correction factor for 
three-dimensional bodies which is independent of the shape of 
the body. It is also generally agreed that the effects of com 
pressibility on the surface pressures of a body in a subsonic 
stream become less pronounced as the body is made more slender 
This latter behavior has been illustrated by examining the maxi 
mum surface velocity in the field of flow about ellipsoids while 
varying thickness ratios, aspect ratios, and free-stream sub 
sonic Mach Number. 

The statement is made by Sauer® that, to a first approxima 
tion, the surface pressures over a slender axially symmetric body 
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of revolution situated in a compressible stream remain unchanged 
from the surface pressures experienced in an incompressible 
stream at the same velocity. This statement is supported by 
in analysis employing appropriate source-sink distributions to 
approximate the presence of a body in the stream. 

The mechanisms of this particular analysis remain somewhat 
obscure. The purpose of this note is to present a clearer illus- 
tration of the statement that to a first approximation there is no 
effect of compressibility on the surface pressures of a slender 
axially symmetric body in a subsonic stream. 

The partial differential equation of continuity for the pertur- 
bation velocity potential in an incompressible stream in cylin- 
drical coordinates is, for axially symmetric flow: 


0°’ 
Ox’? 


0°’ 
or’ 


10’ 


r’or’ 


= () 1) 


where ¢’ = antigradient of perturbation velocities. A general 


solution of Eq. (1) is 


—=—= = 


¢ = — 
V(x’ — 2? +7" 


4x —c/2 


where f(£) is a function representing the strength distribution of 
sources and sinks distributed continuously along the x-axis be- 
2 and the end 


tween the vertex of any given body at x’ = —c, 
of the body at x’ = ¢/2. 
Prandtl and Glauert are credited with introducing the follow- 


ing affine transformation relating the incompressible space x’, 


r’,w’ with a compressible space x, 7, w 


xr = Xx 
on ft 
r ='*/Br = 
p (3) 
Ww Ww 
@ = Ko’ 


where A is any multiplying factor. 

Employing this transformation in conjunction with Eq. (1 
there results the linearized partial differential equation of con 
tinuity for the perturbation velocity potential in a compressible 


stream for axially symmetric flow. 


0" O*¢ ] re) 
32 bs af. : 4. ? = () (4) 
oo“ oF ror 
where B = Y1 — JM? and JJ = free-stream Mach Number. 


Corresponding to Eq. (2), a solution to Eq. (4) using transforma- 


tion (3) is 


— (K/4r) f(g 
Eq. (5) for slender bodies is subject to the boundary condition 
that the flow at the body surface must be tangent to the body 


surface 
or/(U + dr) = (¢7/U)y 0 = tan 0 (6) 
where 6 = angle between tangent to surface of body and x-axis 
and U = free-stream velocity. 
If the given body is to retain its shape no matter what the 
variation of free-stream Mach Number, then 
tan 06 = (¢-/U), —.9 = constant (7 
Returning to Eq. (5) 
K p?r *c/2 f(dé 
or = 8) 
4a _ yo l(x — &)? + Br2]*/2 


This integral is to be evaluated as r approaches a small quantity 
hereafter designated as r — «. This can be done in a more lucid 


manner than the process employed by Sauer.® 


K c/2 


irr otl(x — &) 


f(&) (dé/Br) 
Br\2 + 1}° 


(9) 
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let 
x-—é vx + (c/2 v — (c/2) 
— = ft, = a, —- =) 
Br Br Br 
then 
- *) 
K E P dt 
aS ye — Grt} — x 10 
4rr (t* + L) 7? 
Ja 
Eq. (10) is evaluated now as r — e by taking advantage of the 
mean value theorem. Asr—e,a— ©, b ~ —o, and f(x - 
Brt) — f(x,), where x, is some value of x between x = c/2 and 
e =c/2 
— Kf (x1) wer dt Kf (x1) 
oe er aie 1 
’ trr ([é? + 1] 2ar 
Eq. (11), with eq. (7), yields 
Kf(x:)/2ar = U tan @ 12 


Eq. (12) illustrates that no compensating factor A containing 
Mach Number is needed in the equation to preserve the required 
independence of Mach Number. 

Within the rigor of the linearization implied in Eq. (4), the 
expression for the ratio of the compressible stream pressure co- 
efficient to the incompressible stream pressure coefficient at any 
subsonic free stream velocity U is 


(Cy,)y toes ,= u/u’ 13 


re PO’r re€ 


From Eqs. (2) and (5) for slender bodies with K = 1, 


1 ese + f(t)dé 
P x')r’ +e! = (PxIr +e = ¢)2 ' 
Riis a (X= &)* 
T -/9 
and therefore, 
iC. = (C ’ , 15 
(Com), 4 ¢ Wet oe 


The practical importance of this knowledge lies in the con- 
sideration of the effects of thickness ratio or aspect ratio on the 
Simply stated, the more slender 


compressibility phenomena. 
effects of compressib ility 


the body, the less apparent are the 
Thus the airflow over fuselages will be less influenced by th 
consideration of the compressibility effects than will the airflow 


over wing and tail surfaces with higher aspect ratio. 
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Structural Design for Minimum Weight 


FL. Cox 

National Physical Laboratory, Teddington, Middlesex, England 

April 2, 1949 

— READ with great interest Mr. Shanley’s paper (in th 
March, 1949, issue of the JOURNAL OF THE AERONAUTICAL 

SCIENCES) on the ‘Principles of Structural Design for Minimum 

Weight.””. The conclusions regarding the most economical de- 
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sign for tubes in bending are particularly interesting. It seems 
that, if the allowable stress decreases with ¢/D when that ratio 
tends to zero and if it tends to an upper limit when ¢/D is mod- 
erately large, the optimum working stress and the optimum 
value of t/D depend only on the material and the form of con- 
struction and not at all on the value of the bending moment. 
Moreover, it appears that the optimum stress cannot differ 
much from about half the upper limiting value.. Broad con- 
clusions of this nature are obviously of great value to designers, 
so that it would be worth while to extend the analysis to further 
representative types of construction 

A point that is insufficiently emphasized in Mr. Shanley’s ref 
erence 2 is that the frame spacing in itself does not seriously af- 
fect the structure weight of a stiffened panel—that is, the differ- 
ence in weight between two optimum structures, one with a few 
fairly heavy frames (ribs) and the other with more numerous 
lighter ones, is not marked, If this conclusion may bé applicable 
iso to stiffened tubes in bending, and it seems that it should be, 
the effect of frame spacing in itself may be less than appears from 
Fig. 21. 

It is gratifying to note that the nomenclature and notation 
for analysis of problems in design are becoming standardized. 
For my part I shall be glad to adopt “structural index’’ in prefer- 
ence to ‘structure loading coefficient’’; but I cannot refrain from 
pointing out a lapse by Mr. Shanley himself at the bottom of the 
second column on page 138! 

In conclusion, I must express my entire concurrence in all Mr. 


Shanley’s conclusions and recommendations 


+ 


FORUM 145 


Effect of Rotary Inertia on Higher Modes of 
Vibration 


Henry E. Fettis 
Air Materiel Command, Wright-Patterson Air Force Base. Dayton 
hio 


March 24, 1949 


Te PROBLEM OF THE EFFECT of rotary inertia on the trans 

verse bending vibration modes of aircraft wings is one that 
is frequently neglected when analyses are made. 
reason is the lack of a convenient method for accounting for this 
additional effect. In a recent note (JOURNAL OF THE AERONAU- 
TICAL SCIENCES, Vol. 15, No. 7, p. 425) Professor Scanlan has 
presented a simple and straightforward method of handling this 
In this investigation, as well as in others where the 


Perhaps one 


problem. 
question of rotary inertia is mentioned, the results are confined 
to the effect on the primary mode only, for which the effect on 
frequency and mode shape is negligible even when the inertia is 
extremely large. The reader is thus led to the erroneous con 
clusion that the additional work involved is not warranted in 
view of the slight change that results. However, when the effect 
is included in the second and higher modes, it is found that a sur- 
prisingly large change occurs both in the frequency and in the 
mode shape. To illustrate this, Professor Scanlan’s example is 
employed, and the second modes are computed for varying 


amounts of inertia. 


From the theory of matrices, it is known that the primary mode of the matrix may be suppressed by subtracting from the original 


natrix the matrix 


OF ® 1) ” 

br(2, 4) Nile 

bp(3, 4 ms 
br(4, 4) ms 
5r(4, 4) 


where [ y;')] is the first mode 


By carrying out this reduction process and iterating the second mode 


the following results are obtained: 


I, (y," = 1) 


matrices for various values of the tip weight moment of inertia, 


27,000 Ib.in.* 36,000 Ib.in.? 


] 0 3,000 Ib.in.? 

yn? 2.11 48 0.531 0.464 

yo(? 3.82 82 1.225 1.096 
y,(? 1.00 00 1.00 1.00 

y,‘? —4.98 —3.53 — (0.632 —0.304 

yal? —0.275 —0.111 —(.091 

we 3,310 ¢.p.m. 3,047 ¢.p.m 1,961 ¢.p.m. 1,791 ¢.p.m 
Per cent change 10 40) 46 

The reduction in frequency is now seen to be nearly 50 percent is positive, while for J = 36,000 the tip amplitude is negative, 


for J = 36,000 Ib.in.,? and the effect on the mode shape is also 
large. The deflections are normalized to unity at station 3 rather 
than at the tip, since the tip amplitude approaches zero as J is 
increased, It is, in fact, possible to obtain zero deflection at the 
tip without the presence of an infinite mass, a result that could 
not occur if rotary inertia were neglected. This is brought out 
by computing the second mode for J = 54,000 Ib.in.,? which 
vields the following shape: 


y% = 0.379 
xy = 0.951 
y3 = 1.000 
yw = 0.097 
ys? = —0.066 


w: = 1,575 cycles per min.). In this case all the tip amplitude 


indicating that for some intermediate value of J the tip ampli- 


tude will be zero. 


The Rolling and Yawing Moments Due to Yaw 
of Unswept Wings 


R. F. Anderson 


Northrop Aircraft, Inc., Hawthorne, Calif. 
March 24, 1949 


| is MAKING ESTIMATES of the rolling and yawing moments due 
to yaw of wings with sweep, the increment due to sweep is 
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usually added to the value for the unswept wing. For an un- 


swept wing at an angle of yaw, there is an unsymmetrical lift 
distribution.!. Rolling and yawing moments result, and formu 
las for their values have been obtained by Weissinger.! How 


ever, reference 1 contains errors, and, although corrections were 


made in reference 2, it is not readily available in clearly printed 
form. Certain errors were also made in American reports. It 
is the purpose of this letter to give the correct formulas and to 
make comparisons with some experimental data 


SYMBOLS 


a angle of attack measurement in the plane of symmetry, rad 

¥ = angle of yaw, rad. (see Fig. 1 for positive direction) 

A aspect ratio 

er root chord 

ct tip chord 

r taper ratio ct/er 

K a constant depending on the spanwise flow component (tip shape 

Cr. = lift coefficient 

Cl rolling moment/qSb, stability axes (see Fig. 1) 

C! dC1/dy 

Cl, the rolling moment coefficient at zero lift 

Cm the pitching moment coefficient at zero lift (about the stability 
axis) 

Cy yawing moment/qSb 

Cny WCn/d¥ 

Cng the yawing moment coefficient at zero lift 

m the wing lift curve slope in radians (at zero yaw) 

m the airfoil section lift curve slope, rad. 


THE ROLLING MOMENT 


In carrying out the integrations to find the rolling moment, 
errors were made on page 51 of reference 1. <A value of a con 
stant of 0.29 was obtained instead of 0.15. The constant is given 
correctly in references 2 and 3, and the formula for the rolling 


moment derivative is then 


Cy Kyl + 0.15 [(G-/en) — i —_ : 
f 


oC, A) (cr/cr) + 1 
This is for a wing with linear taper, and A is a constant 
that depends on tip shape (spanwise flow component If 
the taper ratio is now defined as \ = c;/c,, then Eq. (1) may be 
written: 


Cy K/085 + 0.15 


i — 0.05 2 


oC] A 1+hA 
The sign and symbol conventions are shown in Fig. 1 
For an elliptical wing,! 
OCy kK In(4r7A) — 2 
= ().27 = 3 
OC; A 1+ 4 


Weissinger states that the value of A may vary between L| 
and 1.5. For rectangular wings of aspect ratio 5, a value of 1.2 
gave reasonable agreement with tests? (with the yawing moment 
also considered). 

A comparison of the theory and a few tests is shown in Fig. 2 
In the formula, K was taken equal to 1.0 for convenience in com 
puting. The lines are drawn dashed beyond aspects ratios + 
and 12, since the formula is valid only between those limits 
The test data are for the initial straight portion of the C; curves 
i.e., from C, = 0 to about 0.5. The theory and tests are in rea 
sonable agreement if it is considered that for swept wings the 


discrepancies are small compared to the effect of sweep 


THE YAWING MOMENT 


Weissinger’s formula for the yawing moment due to yaw con 
tains a constant 0.02 in reference 1, which is corrected to 0.084 
in references 2 and 3. Conversion from the half span to the 
span gives 0.042, and the formula for the yawing moment deriva 


tive may be written: 
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Pt. No A.R. Ci/ Cy Section® 

1 5.02 1.0 23012 

2 5 1.0 23012 

3 6 1.0 0012 

4 3 1.0 0012 

5 4.13 1.0 23012 

6 5 0.6 23012 

7 5 0.5 0012 

8 5 0.2 23012 

Qg 4.62 0.5 0015 to 23009 


'N.A.C.A. Reports. 
2 Rolling moment, Fig. 16; yawing moment, Fig. 17. 
1 N.A.C.A. 


Wn ; N+ 0.15 


— 0.05 + 0.042 
wc? A 1+ A) ) x 


AVL + V1/a) | (4) 
A+2 m 

This is for a wing with linear taper and no flaps or camber The 

value of m may be taken as the lift curve slope of the unyawed 


wing (in radians). In Eq. (4) the substitution 


aC, = C,?/m 


has been made. 
For the elliptical wing the formula is: 


OCny BP 4 A — 2\1In(4rA) — 2 }1 " 
= —().27 + Z (2) 
oC,” A A+2 A+4 m 
where the sign after A/A has been changed from — to + from 


reference 1 to reference 2. 

A comparison of a few tests with the theory as given by Eq. 
1) is shown in Fig. 3. A value of K = 1.0 was again used for 
convenience. Values of m were obtained from reference + using 
») 


Mo = Z 

In Fig. 3 the agreement between theory and tests is not so 
good as for the rolling moment. This is unfortunate in view of 
the fact that the yawing moment due to yaw of the unswept wing 
is large compared to the increment due to sweep. 

The calculated yawing moment (Fig. 3) does not vary much 
with taper ratio for aspect ratios—say, between 5 and 8. This 
suggests the use of a constant factor for making rough estimates. 
For example, the relation 0C\/ny/0C.2 = —0.125/A could be 
used. This line is drawn on Fig. 3. 

A similar relation due to Hoerner® was given on page 29 of refer- 
ence 6. Even though Hoerner’s formula agrees with some tests, 
its derivation was erroneous. 

The foregoing formulas apply strictly to wings without camber 
For wings with 
The rolling 


or flaps—i.e., to wings for which C,,. = 0. 
camber or flaps, zero lift terms should be added. 


moment increment?:*:7 is given approximately by 
OC, /O¥ ” Ce A (7) 


Fair experimental verification was obtained in reference 2, Fig. 5. 
However, agreement was not so good in reference 3. 

The value of 0C,,./dy is not given explicitly but would appear. 
to be related to the distribution of profile drag along the span. 
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TABLE OF Test DaTA PLOTTED ON Fics. 2 AND 3 


Tip Shape 


Reference 





Remarks 


Revolved T.N. 1046, Fig. 38! No tares 
Revolved Jacobs, U.M. 2069 Table 2, page 16? 
Square RM L7D23' Fig. 36 

Square RM L7D23' Fig. 36 
Revolved T.N. 1046, Fig. 38! No tares 
Revolved Jacobs, U.M. 2069 Table 2? 
Revolved Weissinger, U.M. 1392 Fig. 7 

Revolved Jacobs, U.M. 2069 Table 2? 
Revolved McCormack, RM Fig. 22 


A6K15' 


Elevators on Horizontal Tail Surfaces, N.A.C.A. T.N. No. 1175, February 
1947. 

6 Hoerner, S., 
No. 906. 

6 Toll, T. A., and Queijo, M. J., Approximate Relations and Charts for 
Low-Speed Stability Derivatives of Swept Wings, N.A.C.A. T.N. No. 1581. 

7 Multhopp, H., Die Anwendung der Tragflugeltheorie auf Fragen der 
Flugmechanik, L.G.L. Bericht S.2, Preisausschreiben, 1938/39, 


Forces and Moments on a Yawed Airfoil, N.A.C.A. T.M 


Importance of Intermolecular Forces at Low 
Specific Volumes in the Steady, One-Dimen- 
sional Flow of a Gas 


Lester Lees 

Associate Professor of Aeronautical Engineering, Princeton Uni- 
versity, Princeton, N.J. 

April 4, 1949 


| rHE READERS’ Forum for March, 1949 (pages 190-191), 

Prof. J. D. Schetzer calls attention to the effect of variations 
in the specific heats with temperature on the gas dynamic relations 
for one-dimensional flow. However, in presenting the results of 
the calculations, no mention is made of the effects of intermolecu 
lar forces, as expressed in the dependence of internal energy on 
volume, or the dependence of enthalpy on pressure. Particularly 
at low specific volumes, these effects can easily be large enough 
so that under certain conditions the ‘‘corrections’’ to the perfect 
gas relations are of opposite sign to those predicted by Professor 
Schetzer’s calculations. For example, in an isentropic expan- 
sion the enthalpy change is given by 


dH = c,dT + (0H/d,)rdp = Vdp 


The relation between the differential temperature and préssure 


variations is 


TeV/T), — y 
aT = } a + Lap 


where the quantity in brackets in the Joule-Thomson coefficient. 
Thus, whether the pressure change for a given temperature 
change is smaller or larger than that which would be calculated 
from perfect-gas relations depends upon the sign of the Joule 
Thomson coefficient, and its magnitude relative to V/cp. Simi 
larly, the differential change in volume depends upon the mag 
nitude of the quantity 


cv/{1 + (1/p)(QU/OV)z] 


as compared with the value of cy‘”? for a perfect gas. 
In order to investigate the effects of intermolecular forces, we 


have calculated the isentropic expansion of air from rest for vari 
ous initial conditions, assuming that the gas obeys the Beattie 
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Bridgman equation of state.* In one case the initial stagnation 


temperature is 1,500°R. and the initial pressure 495 Ibs. per sq.in. 
5.0, we find p/po = 0.00187, as compared 


ibsolute; at M = 
= In this 


with p/pp = 0.00189 for a perfect gas with y = 1.4. 
case the “‘correction factor” is 0.989 rather than 0.830, as given by 
Professor Schetzer for M = 5.0, local T = 520°R. The differ- 
ence is explained by the fact that Professor Schetzer’s calcula- 


tions are based on the tables of Keenan and Kaye, which were 


prepared for a pressure of 1 lb. per sq.in. atmosphere and which 
neglect the dependence of internal energy on volume and of en- 


* A report is in preparation 
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thalpy on pressure. At high temperatures cy‘~) and cp age 
both larger than the perfect-gas values; consequently, if one 
neglects the Joule-Thomson effect and puts (OU /0V)7 = 0, the 
volume and pressure changes for a given temperature change 
are larger than that predicted by the perfect-gas relations, and 
the “correction factors” for an isentropic expansion are less than 
=vidently, such a procedure is justified only if the pres- 
One may con- 


unity. 
sures are low—i.e., specific volumes are large. 
clude that no universal correction factors independent of spe- 
cific volume exist (as Professor Schetzer is undoubtedly aware), 
and the one-dimensional flow calculations must be carried out 
for a series of representative values of pressure and temperature, 





Relations Between Viscous and Structural Damping Coefficients 


(Continued from page 410) 


case of Eq. (20), however, the ellipse remains the same 
regardless of changes in frequency. It is ‘a well-known 
fact that in any given structure under static test, as 
the load is slowly cycled through its range of values and 
load-deflection diagrams plotted covering the entire 
cycle, a definite hysteresis loop is obtained which bears 
considerable similarity to the ellipse described by Eqs. 
(20). This hysteresis loop shows no evidence of re- 
ducing to zero as the rapidity of loading is decreased to 
This agrees at least qualita- 
There is also evidence to indi- 


extremely low values. 
tively with Eqs. (20). 
cate that the size of the hysteresis loop does depend on 
the rapidity of loading. However, it is outside the 
scope of the present paper to discuss these deviations 
of an actual structure from the assumed damping 
laws. 

The damping in a structure is frequently defined by 
the ratio of the area of the ellipse given by Eqs. (21) 
and (22) to the strain energy of the structure at its ex- 
treme displacement. The strain energy is (1/2)kx»’, 
and the area ratios are: 


mewxy?/[(1/2)kxo?] = 2rcw/k 


mkgxo?/[(1/2)kxo?] = 21g 


If the damping in a structure does follow the mech- 
anism implied by the assumption of the complex 
stiffness, then the equivalent viscous damping coeffi- 
cient is a linear function of the forcing frequency and 
must increase at the same rate that the frequency de- 
creases. 

Eqs. (23) and (24) define a ‘“‘per cent energy loss due 
to damping,” in that of the total strain energy present 
at an extreme displacement the proportion given by 
Eqs. (23) or (24) is dissipated in completing a cycle 
of vibration. This per cent loss can be a sizable figure 
even for relatively small amounts of damping. For 
example, for a value of g = 0.10 (corresponding to A = 
0.05) 27g becomes 0.628, or nearly 63 per cent of the 


strain energy present at the beginning of a cycle is 
dissipated in the course of that cycle. 


REFERENCES 


1 Theodorsen, T., and Garrick, I. E., Mechanism of Flutter, 
N.A.C.A. Report No. TR685, 1940. 

2 Coleman, R. P., Damping Formulas and Experimental Values 
of Damping in Flutter Models, N.A.C.A. T.N. No. TN751, Feb- 
ruary, 1940. 

3 Den Hartog, J. P., Mechanical Vibrations, 3rd Edition, p. 
53; McGraw-Hill Book Company, Inc., New York, 1947 








are 
if one 
O, the 
shange 
s, and 
s than 
> pres- 
y con- 
of spe- 
ware), 
~d out 


rature 


nech- 
mplex 
-oe ffi- 
y and 


figure 

For 
oA = 
f the 


rcle is 


Flutter, 


| Values 
1 ‘ Feb- 


tion, p 





